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This dissertation is concerned with the dynamic analysis of three- 
dimensional elastic beams which experience large rotational and large de- 
formational motions. To this end, the beam motion is modeled using an 
inertial reference for the translational displacements and a body-fixed refer- 
ence for the rotational quantities. Finite strain rod theories are then defined 
in conjunction with the beam kinematic description which account for the 
effects of stretching, bending, torsion, and transverse shear deformations. 
A convected coordinate representation of the Cauchy stress tensor and a 
conjugate strain definition is introduced to model the beam deformation. 
Due to the inertial reference of the beam kinematics and the convected ref- 
erence of the beam stresses, the present formulation is easily interfaced with 
general multibody dynamics methodologies as well as software modules for 
active control simulations. 

The numerical treatment of the beam formulation is considered in 
detail. A procedure to compute the beam internal force is derived from the 
continuum formulation. The procedure is proven to be invariant to arbi- 
trary rigid motions of the beam while accurately modeling the beam strain. 
To treat the beam dynamics, a two-stage modification of the central differ- 
ence algorithm is presented to integrate the translational coordinates and 
the angular velocity vector. The angular orientation is then obtained from 
the application of an implicit integration algorithm to the Euler parame- 
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ter /angular velocity kinematical relation. The combined developments of 
the objective internal force computation with the dynamic solution proce- 
dures result in the computational preservation of total energy for undamped 
systems. 

The present methodology is also extended to model the dynamics 
of deployment/retrieval of the flexible members. A moving spatial grid 
corresponding to the configuration of a deployed rigid beam is employed as 
a reference for the dynamic variables. A transient integration scheme which 
accurately accounts for the deforming spatial grid is derived from a space- 
time finite element discretization of a Hamiltonian variational statement. 
The computational results of this general deforming finite element beam 
formulation are compared to reported results for a planar inverse-spaghetti 
problem. 
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CHAPTER I 


INTRODUCTION 

The simulation of flexible multibody systems is becoming an increas- 
ingly important tool for the design and operation of many engineering ap- 
plications. The interest in modeling the dynamics of multibody systems has 
emerged in the distinct fields of space dynamics, mechanisms, and robotics. 
Typical examples of these systems include elastic linkages, high precision 
machine dynamics and robot manipulator arms, aircraft propellers, heli- 
copter or turbine rotor blades, and flexible satellites and other types of 
deployable space structures. The articulated structures are thus comprised 
of flexible components which undergo large relative displacements and rota- 
tions in order to carry out the intended operations. To perform the desired 
kinematic motions, various types of mechanical joints are introduced to con- 
strain the relative motion between the various components. 

New technology needs of both the space and robotics industries have 
increased the demand for accurate numerical simulations of the performance 
of multibody systems. The design trend of newly developed mechanisms is 
toward the use of very lightweight structural components. Likewise, equip- 
ment performance requirements are being emphasized which dictate the high 
speed operation and a greater positioning accuracy of these highly flexible 
components. Under these circumstances, a significant coupling can be ex- 
perienced between the gross rigid body motion and the elastic vibrations 
of the mechanism. To accurately simulate this phenomenon and study the 
effect of component flexibility on the overall system performance, a realistic 
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mathematical model of the flexible components that can readily be incorpo- 
rated into a general multibody dynamics methodology is necessary. To this 
end, this dissertation addresses the computational analysis of the dynamics 
of a flexible beam undergoing arbitrary spatial motions and experiencing 
large elastic deformations. 

The development of accurate methods to model the geometric non- 
linearities of structural components has been the subject of numerous inves- 
tigations. One approach considers the small elastic deformations of compo- 
nents which undergo large overall rigid body motions. The two effects are 
modeled separately by introducing a floating reference frame which follows 
some overall mean rigid body motion of the beam; the elastic deformation of 
the beam is then described relative to this moving reference frame [1-13]. In 
this manner, computer codes for analysis of multi-rigid body systems were 
extended to include structural flexibility by superposing existing linear de- 
formation descriptions onto the rigid motions of the floating reference frame 
[11-13], 

The resulting equations of motion of the above approach are in 
terms of a set of reference coordinates representing the motion of the float- 
ing reference frame and a set of relative elastic coordinates representing the 
deformation. For structures consisting of a rigid main body to which flexi- 
ble appendages are attached, the floating reference frame coincides with the 
rotation of the main body [1-2]. For arbitrary configurations in which the 
choice is not obvious, the floating reference frame is constructed to follow 
some mean rigid motion of the flexible body such that the relative defor- 
mation is minimized [3-6]. To determine this frame, constraint conditions 
must be introduced to offset the additional unknown variables of the mov- 
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ing reference; these constraints can be satisfied through a proper choice of 
deformation modes [3-4]. The construction of this mean axis system and 
the appropriate deformation modes has been presented within the context 
of general finite element structural representations [7]. For multibody dy- 
namics applications, the research effort has been toward the construction of 
elastic coordinates which accurately represent the deformation of these joint- 
dominated constrained systems [8-11]. The proper selection of the floating 
reference frame and the consistent set of elastic coordinates is crucial to 
the success of this approach as the component can deform only as dictated 
by the selected modes. The accurate prediction of a proper representation 
remains a difficult challenge in the modeling of flexible multibody systems. 

The initial floating frame approach is limited by an inherent assump- 
tion of linear deformation theory. However, the use of nonlinear theories 
becomes necessary when current space and robotics industry applications 
are considered due to the emphasis toward lightweight and highly flexible 
components. Another instance mandating the use of nonlinear deformation 
theories is the high speed rotation of a flexible component. In this case, the 
rotation gives rise to centrifugal forces which affects the bending stiffness in 
a manner not predicted by the linear deformation theories [14-15]. For these 
purposes, the initial approach was extended to model nonlinear effects by 
including nonlinear strain measures [16-17]. 

An alternative method based on converted coordinate systems has 
been developed to provide the capability of modeling large deformations 
within the context of overall rigid body motions. Instead of modeling the 
entire structural component with a single floating reference frame and an 
appropriate set of deformational coordinates, this method employs finite 
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element based deformational theories by introducing an element-attached 
convected reference frame [18-23]. The convected reference frame is de- 
fined to represent the overall rigid motion of the individual finite element. 
Standard finite element methodologies are then incorporated to describe 
the deformations with respect to the convected coordinate system. This 
method was employed to model large deformations of beam elements using 
updated or total Lagrangian formulations of large deformation continuum 
mechanics [24-27]. A similar concept, termed the element-independent co- 
rotational formulation, was introduced as a method for upgrading existing 
finite elements to allow for large rotations [28-29]. In this method, rigid 
body motions are extracted from the total element displacements prior to 
computing element deformations. Another approach partitioned the compo- 
nent into substructures such that linear elasticity theory referred to a local 
frame was adequate to capture nonlinear effects [30]. 

Recently, a different approach has been adopted to describe the dy- 
namics of a flexible beam which departs from the use of floating or convected 
reference frames [31-36]. The approach introduces finite-deformation rod 
theories from the outset such that the effects of both finite rotation kine- 
matics and large deformations are taken into account [37-43]. The beam 
kinematics are described with respect to the inertial reference frame such 
that the motion due to rigid rotations of the beam is not distinguished from 
that due to deformations. As such, the introduction of a moving frame as 
a reference for elastic deformations is unnecessary. The advantage to this 
is that a natural representation for dynamic systems results such that the 
beam inertia is identical in form to that of rigid body dynamics. In addition, 
large deformations are accurately represented as these formulations incorpo- 
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rate strain measures modeling the combined effects of stretching, transverse 
shear, torsion, and bending deformations. The resulting structure of the 
equations of motion is simpler due to the use of the inertial frame kinematic 
description as opposed to the floating frame or convected frame descriptions. 
These latter approaches result in a complex coupling of the moving refer- 
ence coordinates and the elastic deformation coordinates within the inertia 
terms, whereas the former approach results in a greatly simplified inertia 
operator. In exchange for this simplicity, a nonlinear internal force expres- 
sion which is invariant to rigid body motions must be derived, as existing 
linear deformation descriptions referenced with respect to a fixed coordinate 
system become invalid for systems exhibiting large rotational motions. 

The development of a computational procedure capable of real-time 
simulations of space and robotic applications requires further research in- 
volving both the model formulation and computational solution procedures. 
To this end, the present work is focused on achieving effective and accurate 
computational methods for the simulation of multibody systems in which 
flexible components may be conceptualized as spatial beams. For a compu- 
tationally effective formulation, an automatic derivation of the equations of 
motion is necessary which includes a realistic modeling of the flexible mem- 
bers and a streamlined incorporation of system constraints. Likewise, robust 
and efficient time integration procedures which exploit inherent character- 
istics of the formulation axe required. Finally, it is desirable to be able to 
incorporate additional analysis capabilities such as the deployment and/or 
retrieval of the flexible member, active control and state estimation, ther- 
mal and environmental effects, and other specialized fields into the present 
formulation. As each of these computational elements are best formulated 
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by specialists of the various fields, it is advantageous to conduct a modular 
solution approach toward the software development rather than embedding 
several analysis capabilities into a single monolithic program. 

To this end, separate software modules have been developed which 
are easily interfaced to form an attractive multibody dynamics simulation 
package. First, a computer-oriented formulation of flexible spatial beams has 
been developed to constitute an integral kernel of the multibody dynamics 
methodology. The flexible beam model accounts for both finite rotations 
and large deformations, and can be incorporated into general multibody 
dynamics systems in a straightforward manner. The beam kinematics are 
referenced directly to the inertial frame. As such, the equations of motion 
for an arbitrary configuration of flexible beams and rigid bodies can auto- 
matically be generated in terms of an identical set of physical coordinates. 
The structure of the equations, aside from an expression for the nonlinear 
internal force, is identical for both the rigid and flexible components. Nu- 
merical strategies developed for the solution of equations representing spa- 
tial kinematic systems are thus applicable for the entire flexible multibody 
system. Such a unified treatment is not applicable for equations formulated 
within the context of the floating frame approach. In this case, the reference 
and elastic coordinate definitions are of highly different character and thus 
require separate numerical treatment. 

A multibody dynamics solution procedure, originally demonstrated 
on rigid body systems in previous studies [44], has been adopted for the 
present flexible multibody formulation. For multibody dynamics applica- 
tions, time integration algorithms must include a proper treatment of three- 
dimensional finite rotations and also a method to satisfy the kinematic con- 
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straint conditions. In the present work, an attractive partitioned solution 
procedure has been developed such that the generalized coordinates are 
solved in a separate module from the constraint forces. The generalized 
coordinate solution procedure is based on the application of an improved 
variation of the explicit central difference algorithm to the translational co- 
ordinates and the angular velocity. To update the configuration orientation 
from the angular velocity, a separate numerical procedure based on the Eu- 
ler parameter representation of finite rotations is introduced. This overall 
solution of the spatial dynamics is successfully interfaced with a separate 
module which enforces the multibody system constraints. The constraint 
force solver implicitly integrates a stabilized companion differential equa- 
tion for Lagrange multipliers. 

The application of these multibody dynamics solution procedures 
to systems including spatial beam components relies on an accurate com- 
putation of the beam internal force. As the degrees of freedom implicitly 
contain information of both rigid motion and strain deformation, the inter- 
nal force computation depends on a judicious procedure which filters out 
the rigid body motions embedded within these variables. For this purpose, 
an objective strain increment/stress update procedure has been developed 
which remains invariant to arbitrary rigid body motions occurring in finite 
increments of time. The combined developments of this internal force com- 
putation with the multibody dynamics solution procedures results in a com- 
putational preservation of total energy for undamped systems. This distinct 
feature of the present work will be demonstrated with several examples. 

Simulations of active control/ vibration suppression or controlled 
slewing maneuvers can easily be performed using the present beam for- 
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mulation. This is due to the use of a convected coordinate representation of 
the Cauchy stress tensor to describe the beam deformation. The physically 
appealing description of stress as related to the surface of the deformed con- 
figuration has been recast into a coordinate system moving with the beam 
in a manner which provides conceptual and computational simplifications 
of the deformation representation. As such, the actual strains measured by 
sensors located and operating on the deformed structure correspond to those 
with which the present multibody dynamics formulation is based. Control 
software modules may thus be interfaced with the present multibody dy- 
namic solution modules in a straightforward manner. 

The present methodology is also extended to include the dynamics 
of deployment and/or retrieval of the flexible members. Hamilton’s law of 
varying action is used to formulate the equations of motion for this three- 
dimensional moving boundary value problem. A moving node reference for 
the beam dynamics is employed within the present formulation to account 
for the changing spatial volume. The reference, which previously was fixed, 
now corresponds to the configuration of a deployed beam as if it were rigid. 
A transient integration scheme is derived from a space-time discretization of 
the Hamiltonian formulation, effectively accounting for the changing spatial 
reference. The methodology is successfully interfaced with the internal force 
computational procedure, thus retaining the large rotation/large deforma- 
tion modeling capabilities of the present work. The formulation and the 
computational procedure are then specialized to a planar inverse-spaghetti 
problem for illustrative purposes. 

The rest of the thesis is organized as follows. Necessary mathe- 
matical preliminaries are presented in Chapter II in deriving the equations 
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of motion for a flexible spatial beam. The numerical treatment of a dis- 
crete form of the equations is then considered in detail in the remaining 
chapters. A computational procedure for the internal force is derived from 
the continuum formulation in Chapter III. Special care has been taken to 
achieve a computation which remains invariant to arbitrary finite rigid rota- 
tions while accurately representing the strain of the beam. The formulation 
is then extended into the multibody dynamics framework in Chapter IV. 
The solution techniques for multibody systems, including the orientation 
update procedure and the treatment of constraints, are presented in this 
chapter. The combined effort of these techniques with the invariant inter- 
nal force computation results in successful simulations of flexible multibody 
dynamic systems as shown in several examples. The extension to model de- 
ployment/retrieval dynamics of the beam member is presented in Chapter 
V. The methodology developed to simulate the dynamics of the changing 
spatial volume is presented in this chapter. Chapter VI then summarizes 
the major contributions of this work. 




CHAPTER II 


FLEXIBLE BEAM DYNAMICS FORMULATION 
2.1 Introduction 

The formulation of the equations of motion for a flexible spatial 
beam is presented in this chapter. The formulation accounts for the effects 
of both finite rotations and large deformations of the beam component. 
To model these geometric nonlinearities, a relevant description of the spa- 
tial kinematics must be introduced. The equations of motion of the beam 
as governed by the basic principles of continuum mechanics are then de- 
rived from these kinematic definitions. A complete understanding of this 
continuum formulation is necessary prior to the development of effective 
computational techniques for a finite strain beam theory. 

To formulate the equations, the beam kinematics are described with 
reference to the inertial frame as in [31-36]. With this approach, the motion 
due to rigid rotations of the beam is not distinguished from that due to de- 
formations. The beam configuration is defined by a position vector locating 
the centroid of a typical cross-section and an orthogonal matrix designating 
the orientation of this cross-section, both of which are referenced to a fixed 
inertial frame. As a consequence, the beam inertia operator is identical in 
form to that of rigid body dynamics. This procedure fully departs from the 
traditional approach discussed in the preceding chapter in which a floating 
reference frame is introduced to separate the rigid body motions from the 
elastic deformations. 
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As the kinematics are referred to a fixed inertial coordinate system, 
the present formulation relies on the use of strain measures which are invari- 
ant to rigid body motion. Applicable invariant strain measures which model 
stretching, transverse shear, torsion, and bending deformations of rod-type 
structures have been developed in the literature. The classical Kirchoff-Love 
rod theory modeling large bending and torsion deformations was extended 
to include stretching and transverse shear strains in [37]. Other formula- 
tions derive similar rod theories from an appropriate version of the principle 
of virtual work [38-40] and from kinematical considerations including the 
effect of warping which induces torsion-bending coupling [41-43]. The use 
of these rod-type theories within dynamics problems involving laxge spa- 
tial rotations was initiated in [31-33] and further pursued in [34-36]. These 
and other works modeling the finite rotations and deformations of beam 
components employ the Piola-Kirchoff stress representations in which trac- 
tion forces acting on a deformed surface element are referenced back to the 
undeformed configuration [24-27,31-36]. Although these formulations are 
mathematically consistent, the stresses transmitted in the instantaneous 
state are referred to the initial state in a way that is physically artificial 
[45]. Transformations must then be made to relate the stresses back to the 
actual deformed configuration. In addition, the use of the unsymmetric first 
Piola-Kirchoff stress tensor as employed in [31-36] leads to complexities in 
subsequent computations and linearizations. 

A desire for a beam formulation which can easily be interfaced with 
active feedback control schemes has motivated a more physically-based in- 
terpretation of the spatial beam formulations. To achieve real-time software 
simulations, the computed deformation representations must correspond to 
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the actual stress/strain measurements of the sensors located and operat- 
ing on the deformed structure. For this purpose, the present formulation 
employs the Cauchy stress and a corresponding strain tensor, which are 
related to the true deformed beam configuration, to describe the elastic de- 
formations. To implement this physical representation, a convected frame 
which follows the rigid motion of the beam is introduced as the reference for 
the Cauchy stress and corresponding strain components. Thus the present 
formulation combines this convected reference for the stress representation 
with the inertial reference for the beam dynamics. As such, the present 
formulation is easily adapted within general multibody dynamics method- 
ologies; this advantage will be illuminated in Chapter IV . The formulation 
can also be readily extended to model the dynamics of deployment and re- 
trieval of the beam as will be shown in Chapter V. The key to the success of 
the present formulation within these applications is a procedure developed 
for the computation of the internal force. This algorithmic development, 
presented in Chapter III, guarantees that the physical stress representation 
computationally remains invariant to arbitrary rigid-body motions. 

The rest of this chapter is organized as follows. Section 2.2 will 
detail the beam kinematics in which the total motion is referred directly to 
the inertial reference frame. The description is compared to the kinematic 
formalisms of the floating frame approach, thus illuminating advantages to 
the former description. The equations of motion are formulated in Section 
2.3 by specializing the principle of virtual work of a continuum solid to the 
spatial beam kinematic assumptions. The convected frame decomposition 
of the Cauchy stresses is introduced, and virtual strain-displacement rela- 
tionships are identified from the virtual work expression. The convected 
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coordinate stress/strain representation provides conceptual simplifications 
in the derivation and subsequent computations of the internal force. To 
complete the formulation, a constitutive law relating the convected frame 
stress and strain rates is presented in Section 2.4. Consistent linearization 
procedures employed to obtain linearized weak forms of the equations of 
motion axe then presented in Section 2.5. Using a natural continuum based 
approach for the linearization procedure, the present formulation inherently 
based on stresses referred to a deformed configuration results in symmetric 
tangent matrices. A summary of the present formulation is then given in 
Section 2.6. 

2.2 Large Rotation Beam Kinematics 

The inherent difference of the floating frame approach and the in- 
ertial reference approach to describe the dynamics of spatial beams is il- 
luminated as follows. As stated in Chapter I, the floating frame approach 
introduces a reference frame to follow an overall mean rigid-body motion of 
the beam; the elastic deformations of the beam are then described relative 
to this moving reference. In this manner the motion due to the rigid rota- 
tion of the beam is separated from the local deformation of the beam. The 
position vector describing the location of an arbitrary particle point on the 
beam from the inertial origin, as shown below in Figure 2.1, is given by 

r = z T e + ( X + u f ) T f . (2.2.1) 

The following notation is introduced in (2.2.1) to describe these kinematics: 
e = { ei,e 2 ,e 3 } r represents the orthogonal basis vectors of the inertial 
reference frame; f = { fi,f 2 ,f3 } T represents the orthogonal basis vectors 
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Figure 2.1 Spatial Beam Kinematics: Floating Frame Approach 

of floating reference frame; z T = { Z] , z 2 , Z 3 } represents the inertial com- 
ponents of the origin of this moving reference frame; X = { X 2 ,X 2 , A 3 } T 
represents the moving frame components of the undeformed beam position, 
and Uf = { Uf lt Uf 2 ,Uf a } T represents the moving frame components of 

the relative deformational displacement of the beam from the undeformed 
position. The terms z T e and X T f represent the rigid-body motion of the 
beam, while the terms u J f represent the small deformational displacements 
relative to the rigid-body reference frame. 

The 3x3 orthogonal transformation matrix A is introduced to orient 
the f-basis from the e-basis as f = A e. The angular velocity of the moving 



15 


frame is thus given by [47-48] 
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in which the tensor components are referred to the moving basis and the no- 
tation [ ' ] represents skew-symmetric matrices. A conjugate virtual rotation 
tensor is defined analogous to the angular velocity tensor as 


60 = —6AA r 


(2.2.3) 


The time differentiation of (2.2.1) is obtained from the well known formula 
applicable for rotating coordinate systems [47] 


— = — - — 
dt ~ dt dt 


(2.2.4) 


in which u? is the angular velocity vector and the superscripts e and b indicate 
that the derivatives are to be those observed in the inertial and moving 
system of axes respectively. The above relation is expressed in the following 
matrix form as 


± _ d*_ 

dt dt 


(2.2.5) 


which acts on the moving frame components of a given vector. The velocity, 
acceleration, and variation of the position vector (2.2.1) axe thus given as 


r 

« T 

— z e 

+ 

d T f 

+ ( 

X T + u T ) uj f 

(2.2.6a) 

r 

-T 

= z e 

+ 

ii r f 

+ 2 

u T uj f -I- ( X T + u T ) 



rT f 

Wj I 

+ 

{X T 

+ « T 

)uj uj { 

(2.2.66) 

Sr 

= 6z t e 

+ 

6u t 

f + 

( X T -(- u T ) 60 t f . 

(2.2.6c) 


The major drawback of the floating frame approach is immediately seen 
in the above acceleration expression. The terms iif and u / representing 
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the elastic deformations axe coupled with the terms Uf and Uf representing 
the angular velocity and acceleration of the floating reference frame. This 
coupling between the rigid body motion and the elastic deformation requires 
the development of specialized numerical procedures for the solution of the 
equations of motion formulated from this kinematic description. 

The present formulation adopts an inertial reference frame for de- 
scribing the translational motions and a body-fixed reference frame for de- 
scribing the rotational motions. The consequence of this description is that 
the translational and rotational variables embody information due to both 
rigid-body motions and deformations of the beam. The configuration of 
the beam is completely characterized using a position vector locating the 
neutral axis of the beam from the inertial origin and a body-fixed reference 
frame representing the orientation of the cross-section with respect to the 
inertial reference frame. The position vector describing the location of an 
arbitrary particle point on the beam from the inertial origin, as shown below 
in Figure 2.2, is given by 

r = ( X + u ) T e + i T b . (2.2.7) 

In the above equation, the notation e = { ©1,62,63 } represents the or- 

T 

thogonal basis vectors of the inertial reference frame; b = * { b 1? b 2 , b 3 } 
represents the orthogonal basis vectors of the body-fixed reference frame 
attached to the beam cross section; X = { X\, X21X3 } T represents the 
inertial components of the original neutral axis position; u = { , u 2 , u 3 } 

represents the inertial components of the subsequent total translational dis- 
placement of the neutral axis, and i T = { 0,f? 2 >^3 } represents the body- 
fixed components of the distance from the beam neutral-axis to the ma- 
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Figure 2.2 Spatial Beam Kinematics: Nonlinear Continuum Approach 

terial point located on the deformed beam cross-section. Thus the total 
large rotation motions are modeled by a combination of the translational 
displacements u and the moving reference frame b. It is noted that the 
beam cross-section is allowed to rotate with respect to the neutral axis of 
the beam. In this manner, transverse shear deformations are modeled as the 
orientation of the beam cross-section is not necessarily perpendicular to the 
beam neutral axis. Waxping deformation of the cross-section is not taken 
into consideration in the present formulation. 

The orientation of the body-fixed reference frame is expressed with 
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respect to the inertial reference frame as 

b = Re 


( 2 . 2 . 8 ) 


where R is a 3 x 3 orthogonal transformation matrix. The body frame 
components of the angular velocity tensor and the virtual rotation tensor of 
the cross-section orientation are given by 


w = -— R r , 8a = -8 KK t , (2.2.9) 

dt 

respectively. The velocity, acceleration, and variation of (2.2.7) are thus 

* = — e + Fu; T b (2.2.10 a) 

dt dt 

£l = e + e ( — + H T Z T ) b (2.2.106) 

dt 2 dt 2 { dt ' 

Sr = Su T e + e T Sa T b . (2.2.10c) 


It is seen in this acceleration expression, in contrast to (2.2.6b), that the 
translational displacements of the neutral axis are completely decoupled 
from the angular orientations of the cross-section. The acceleration is of 
the same form as in the Euler equations for rigid body motion. This leads 
to an effective partitioned numerical solution procedure which is equally 
applicable to both rigid and flexible components of a multibody system. 
The derivation of the equations of motion are discussed next. The numerical 
techniques for the solution of these equations of motion are discussed in 
Chapter IV. 

2.3 Nonlinear Equations of Motion 

The conditions for the equilibrium of a solid continuum are given 
by Cauchy’s equations of motion. These differential equation of motions are 
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stated as [45-46] 

do^i 

+ >>(f‘ - i >) = 0 

a h = °)i 

to be satisfied through the interior of the continuum V subject to displace- 
ment and traction boundary conditions 

Xi = Xi on S u , ti = (TjiUj on S a 

on the respective displacement and traction surfaces 5 U , S a . The Cartesian 
coordinates x, represent the particle position after the deformation, v, the 
particle velocity, fi the external force per unit mass, and p the mass per unit 
volume. Likewise, a]- represents the Cartesian components of the Cauchy 
stress tensor and ti the stress vector acting on a surface with outward normal 
components n,. 

The weak or variational form of Cauchy’s differential equations 
are employed to deduce a conjugate set of Cauchy stress/virtual strain- 
displacement relations as well as to provide a basis for the displacement 
based finite element method. The variational form of Cauchy’s equations is 
the principle of virtual work given as [45] 

L 6n * dv + X < lx jv = 

(2.3.1) 

The Cartesian coordinates Sri represent a kinematically admissible vir- 
tual displacement. This virtual work expression is tailored to the contin- 
uum beam by incorporating the kinematic relations (2.2.7), (2.2.10b), and 
(2.2.10c) for the components Xi, fi, and Sr,, respectively. For notational 


X 


Sr,f, dV + 


L 


Sr,t t dS 
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convenience and subsequent finite element discretization, the principle of 
virtual work is expressed in the following operator form: 

8F 1 + 8F S = 8F e + 8F t (2.3.2) 

where the inertia operator 8F 1 , internal force operator 8F S , external force 
operator 8F E , and traction operator 8F T are identified from (2.3.1). Ex- 
plicit expressions for the various operators are given in the following Sections 
2.3.1 to 2.3.3. 

2.3.1 Inertia Operator 

The inertia operator defined in (2.3.1) 

8F 1 = [ p 8ri?i dV = f p 8r • r dV (2.3.3) 

Jv Jv 

is tailored to the present problem with the kinematic equations (2.2.10). 
The following simple expression results for 8F 1 if the origin of the body- 
fixed basis is located at the centroid of the cross-section: 

^ • (2 ' 3 ' 4) 
where 

m = J p dA , J = J p dA 

represents the area and inertia properties of the beam cross-section and s 
represents a length parameter to be taken along the beam neutral axis. It is 
seen that the translational inertia is completely decoupled from the rotary 
inertia and is of the same form as the classic Euler equations of rigid body 
dynamics. This is due to the dual choice of the translational displacements 
measured in the inertial basis and the angular velocity measured in the 
body-fixed basis located at the center of mass of the cross-section. 
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2.3.2 Internal Force Operator 

The internal force operator is defined from (2.3.1) as 

iFS = L < 2 - 3 - 5 ) 

identifying as conjugate quantities the virtual displacement gradient and 
the Cauchy stress tensor. An explicit expression for the internal force is 
then derived from the beam kinematics (2.2.10). A set of virtual strain- 
displacement relations that axe invariant to rigid body motions axe deduced 
from this virtual work expression as in [38-40]. To complete the formulation, 
a suitable constitutive relation is chosen to relate the strain rates, which axe 
of identical form as the virtual strain relations, to an appropriate stress 
rate tensor. An objective incremental procedure is then derived from this 
rate-type constitutive law. 

To provide conceptual simplifications in the derivation and subse- 
quent computations, the Cauchy stress tensor and the virtual strain tensor 
axe decomposed with respect to an alternative beam reference frame which 
lies tangent to the deformed neutral axis. The virtual strain tensor contains 
three independent, non-zero components when referenced to this convected 
coordinate system. In addition, the task of stress update is accomplished 
with a much simpler computation when convected frame as opposed to in- 
ertial frame stress components axe considered. When inertial frame compo- 
nents axe employed within problems exhibiting laxge rigid body motions, an 
appropriate nonlinear constitutive relation is required to obtain an objective 
stress update from strain measures. When expressed using convected frame 
stress and strain tensor components, the constitutive law leads to a simple 
additive update procedure. 
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ei 


Figure 2.3 Spatial Beam Kinematics: Convected Reference Frame 

For this purpose, we introduce a convected reference frame which 
lies tangent to the deformed neutral axis. The convected reference frame 
a = { ai,a 2 ,a 3 } T , as shown above in Figure 2.3, is related to the inertial 
reference frame e by a = T e . For present implementation purposes 
within the context of a constant-strain finite element, the convected frame 
will be approximated as a straight line connecting the two element nodes. 
As such, the reference frame is constant on the element level and is simi- 
lar in concept to that introduced in [19-20]. It is noted that this reference 
frame does not coincide with the body frame b attached to the cross-section. 
The relative difference between these two frames is represented by the ro- 
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tation matrix S which models the effects of transverse shear and torsion 
deformations as 

b = S a , R = ST . (2.3.6) 

As will be shown, the interdependence between R, S, and T plays a key 
role in the algorithmic implementation of the present formulation. 

The internal force operator, originally characterized by inertial 
frame components of the Cauchy stress tensor cr^ and conjugate virtual 
displacement gradient, can equivalently be expressed in terms of the con- 
vected frame stress tensor components a fj and a corresponding convected 
virtual displacement gradient. As the rotational tensor T tJ maps the in- 
ertial coordinates x; to the convected coordinates the following tensor 
transformations 


£- = T ti A at s = T ti T,, (2.3.7) 

are incorporated into (2.3.5) to yield 

/ d6r 

T m i dV . (2.3.8) 

The virtual strain tensor Se^k 1S defined from the symmetric portion 
of the transformed deformation gradients as 


6e 


a 

mk 


1 

2 


( T m i 


dSri 

dik 


+ 


T kt 


dSri 

dtm 



(2.3.9) 


The virtual strain definition (2.3.9) vanishes during a rigid body rotation. 
As the virtual displacements of a rigid-body rotation are <5r = £ar, the 
virtual deformation gradient reduces to the skew-symmetric matrix Sot in 
this instance. As such, the symmetric virtual strain definition (2.3.9) is 
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invariant to rigid, body motions and is thus classified as an objective rate- 
type tensor. 

Due to the symmetry of the Cauchy stress tensor, the virtual work 
of the internal force is given as 


6F- 


-L 


femk dV • 


(2.3.10) 


This expression is rewritten in the following vector product format as 

(2.3.11) 

l J 

In the above equation, the notations 


SF i 


r - „ f fe «l 

= / { a (v [ dV ■ 

Jv ( J 


6 = {6,6,6} = {6*7,0 , Oij = ( ),j + Oji 


denote the coordinates of the convected reference frame and the engineering 
shear strain definitions respectively. The remaining convected frame strain 
components 3X6 identically equal to zero due to the orig- 

inal assumptions of the beam kinematics. From (2.3.9), the three non-zero 
virtual strain components are given as 


6e^ 

= T u 

6s^ 

= T 2i 


= T 3i 


dSri 

dSri 

dSri 


+ 

+ 


T u 


T u 


dSri 
d t) 
dSri 

W 


(2.3.12a) 

(2.3.126) 

(2.3.12c) 


The above definitions are rewritten terms of the virtual translations and 
virtual rotations by performing the parametric differentiation of the vari- 
ation of the position vector (2.2.10c) with respect to the convected frame 
coordinates. As the convected reference frame has been constructed to be 
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constant over the element domain where the differentiation is performed, it 
becomes convenient to express the variation of the position vector as 


6r = Su T e + 60 T l J a 

= ( Su T + 60 T i a T)e . (2.3.13) 


In the above equation, the components 60 and t a axe defined as 


60 = S T 6a , l a = S T l = l 


f6 

n 

lc 


(2.3.14) 


to represent convected frame decompositions of the cross-section virtual 
rotation vector and material position coordinates. The variation (2.3.13) is 
differentiated with respect to the neutral-axis coordinate £ to give 


d 6r 

w 


= ( 


+ ?S t L t) 


where the approximation 


dt 

dl a 

dd 


dd 


(2.3.15) 


0 


is employed as the change in the neutral-axis coordinate of the total position 
due to a small rotation of the cross-section about the neutral axis is negligi- 
ble. Likewise, the differentiation of (2.3.13) with respect to the coordinates 
perpendicular to the neutral- axis yields 


d 6 r 

8 t ) 

d 6r 

where and are given by 


60 T i 2 T e 
60 T ijTe 


i2 


0 0 1 

0 0 0 

-10 0 


dja 

dr] 


(2.3.16) 

(2.3.17) 
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*3 


0-10 

1 0 0 

0 0 0 


di a 

dC 


It is noted that in the above derivations the assumptions intrinsic to the 
beam formulation in which the virtual translations 8u and virtual rotations 
< 5/3 are independent of the cross-section variables r/ and ( have been effected. 
Substitutions of the derivatives (2.3.15), (2.3.16), and (2.3.17) into (2.3.12) 
yield the expression 


<5 ^ ij 


= T 


d 8u 


+ t'i <5/3 + i , 


IT 


d8/3 

d( 


(2.3.18) 


where 


1 1 = 


0 0 0 

0 0-1 

0 1 0 


The above equation (2.3.18) is rewritten as 

f 1 

< Sifr, > = 6j + l 6k 


(2.3.19) 


l J 

where £7 represents the membrane and two transverse shear virtual strains 
and 6k represents the torsion and two bending virtual strains. The defini- 
tions 


r „ <9(5 U 

6 *f = T -777- + ^ 


8k = 


di 

d 8/3 
dt 


0 

-8/3 3 
{ 8/3 2 J 


(2.3.20) 

(2.3.21) 


are a convected coordinate representation of the virtual strains derived in 
[39-40]. 

To determine the stress state, it remains to introduce an appropriate 
constitutive law. A rate-type constitutive law, as discussed in Section 2.4, 
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is employed for this purpose. The resulting stress state, as derived in the 
following section, is written in form similar to (2.3.19) as 

f 1 t 

I = <7 7 + F <t k (2.3.22) 

where <7 7 represents the membrane and two transverse shear stresses and 
o K represents the torsion and two bending stresses. From this result, the 
virtual work of the internal force can be written as 


6F S = J { 6j t N y + Sk t M k } di (2.3.23) 

by substituting (2.3.22) and (2.3.19) into (2.3.11) and integrating over the 
area coordinates of a symmetric cross section. In the above equation, iV 7 
represents the axial and transverse shear forces per unit length, and M K 
represents the torsional and bending moments per unit length as 

N 1 = j <7 dA , M k = f l cr dA . (2.3.24) 

Ja Ja 

It is necessary to rewrite (2.3.23) in a manner consistent with the 
inertia operator (2.3.4) as 

SF s = j {Su T fa 3 -} [ B \ T j"’} d( (2.3.25) 

To identify the proper strain-displacement matrix [ B ], a transformation 
of the virtual strains back to the body frame components of the virtual 
rotations is required. To effect the change of the body reference frame of 
the cross-section orientation in space with respect to the constant convected 
reference frame, we invoke the following relations 


d a 6(3 

dt 


c t d a 6a 
T d^Sa 

^ V Ac 


(2.3.26) 


rp 8 a S rp 

*s = —S r (2.3.27) 




-I- k s 6a ) , 
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which axe completely analogous to the time derivatives of rotating coordi- 
nates. The strain operator [ B ] is then identified as 


; b ; 


[ i. S T l 

0 S T ( ks + I ) . 


(2.3.28) 


To complete the internal force definition, it remains to provide a procedure 
for updating cr-y and <r K in order to compute Nf and M K . This is accom- 
plished via the numerical integration of an appropriate rate-type constitutive 


law to be discussed shortly in Section 2.4. 


2.3.3 External Force and Traction Operator 

The external force operator defined in (2.3.1) as 


6F e 


has the final resultant form 



Sr, ft dV 


6F e = j { Su T 6a T } | f f l } d£ (2.3.29) 

where f e represents the inertial components of a force per unit length acting 
on the beam neutral axis and f b represents the body-fixed components of 
a moment per unit length acting on the beam cross-section. The traction 
operator defined as 

SF T = f Sri t { dS (2.3.30) 

Js 

acts on the exterior surfaces of the beam as natural boundary conditions. 


2.4 Constitutive Equations 

To complete the description of the nonlinear continuum problem, 
the derivation of a stress-state from an applicable constitutive equation is 
necessary. The classical elastic constitutive equations applicable within the 
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finite deformation context relates the second Piola-Kirchhoff stress compo- 
nents to a linear combination of the Lagrangian finite strain components 
within a generalized Hooke’s Law [45]. This type of deformation description 
is based on a direct comparison of the the current shape of the continuum 
back to its original reference configuration. Such a formulation thus asserts 
that regardless of the extent of the deformation, the material response of the 
continuum is based solely from the unstressed state without any reference 
to the history of the deformation. A more natural concept is based on the 
idea that an increment of stress is a function of the increment of strain from 
an immediately preceding state [49]. The constitutive law which mathemat- 
ically generalizes this concept relates the instantaneous rate of stress to the 
instantaneous rate of deformation. This concept is adopted in the present 
beam formulation. 

When stress and strain rates are used with a constitutive law, the 
property of objectivity must be taken into consideration. The principle of 
objectivity requires that intrinsic physical properties of a body be inde- 
pendent of the body’s orientation in space. This principle is embodied in 
constitutive theory by requiring that constitutive equations contain only ob- 
jective tensor fields. Tensors defined from the material coordinates of the 
undeformed continuum such as the Piola-Kirchoff stress tensors automati- 
cally possess this property as do time independent tensors. However, the 
time derivative of the Cauchy stress tensor is not objective, and an alter- 
nate rate definition other than the time derivative must be used within a 
constitutive law based on Cauchy stress rates. 

One objective stress rate constructed from the Cauchy stress tensor 
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is the Truesdeil stress rate tensor a k i defined as 

&kl = **tl ~~ a km v l,m ~ °ml v k,m + a kl v m,m > 

and an objective constitutive law incorporating this definition is 


(2.4.1) 


&kl ~ c ijkl £kl 


(2.4.2) 


The notation 

e din , e 1 / dii k diii \ 

U,,m ~ dx m ' £kl ~ 2 dxi dx k 

represents the velocity gradient tensor and the symmetric part of the ve- 
locity gradient respectively, and Cijkl is the spatial transformation of the 
material response tensor Cabcd to be defined shortly. Similar objective 
stress rates, as the Jauman stress rate tensor, may also be employed within 
objective constitutive laws. In the present formulation, the Truesdeil rate 
equation is chosen as it leads to symmetric tangent stiffness matrices within 
the linearization procedure dicussed in Section 2.5. 

For computational convenience, the Truesdeil constitutive equation 
(2.4.2) is rearranged as 

*ij + erf* ~ ™ C ,k a kj = C ljk i i%i (2.4.3) 


where 


Cijkl = c ijkl — SklGij + 2 ( + f>jl&ik + k<?jl + &jk&il ) (2.4.4) 

corresponds to a stress-dependent constitutive tensor and 


- 1 / diik 

w kl = o ( 


dii 


2 v dxi dxk 


L ) 


(2.4.5) 
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corresponds to the skew- symmetric spin tensor. This form (2.4.3) is con- 
ducive for the development of a computational procedure to update a given 
stress state. Any stress rate defined from a skew symmetric tensor as 

<j <w> = a + aCl - Cla , ft = W W T 

is automatically objective. It is easily shown that when this stress rate 
is transformed to the basis defined by the orthogonal rotation tensor W 
characterizing ft as 


a 


<vv> 

w 


W T ( a + aft - fta ) W , 


it is precisely the time derivative of the stress tensor also transformed to the 
W basis as 

°* w> = 5( w3Vw ) = w • 

The transformed stress rate tensor can then be directly integrated as 

tf'vv" 1 = &y/ + f {Cs)w dt (2.4.6) 

Jt n 

where the right hand side of the constitutive equation is also transformed 
to the W basis. The above equation is accurately approximated with a 
midpoint integration rule as [50] 

o’ "v 1 = a w 4" C As w (2.4.7) 

where the increment Aew is strictly due to deformation and calculated at 
time 


This stress update procedure, as applied to the Truesdell rate defi- 
nition, requires that the stress and strain tensors originally referred to the 
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inertial basis be transformed to a basis rotating with the spin tensor. This 
basis must be calculated from its generating differential equation 

W = w W , 

and for this purpose algorithms have been developed which computationally 
remain objective [51-53]. However, this additional complexity is bypassed 
within the present beam formulation. The update procedure (2.4.7) can be 
directly applied to the convected frame tensor components of the Cauchy 
stresses and strains without the need for further transformations. This is 
shown as follows by expressing the the Truesdell rate equation (2.4.3) with 
convected frame tensor components as opposed to inertial frame compo- 
nents. 

To derive the appropriate constitutive law for the present beam 
formulation, the origin of the Truesdell rate equation is first examined. The 
Truesdell rate equation is based on the classical hyper-elastic constitutive 
equations [45] 

( Skl ) = Sab = Cabcd Ecd (2.4.8) 

where Sab is the second Piola-Kirchhoff stress tensor, Ecd is the Lagrange 
strain tensor and Cabcd is a material response tensor derived from a pre- 
scribed elastic strain energy function. The second Piola-IvirchhofT stress 
tensor is related to the Cauchy stress tensor a £ by 

Sab == j %A,i &ij Xbj > 
and the Lagrange strain tensor is defined as 

EaB = \ ( Xi,AX,,B ~ 6 ab ) 


(2.4.9) 
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In the above equations, the components Xj denote the material coordinates 
of the undeformed configuration, x,- the spatial coordinates of the deformed 
configuration, and 


,A — 


dx i 

Wa 


j = det( x itA ) 


X A ,i 


dX A 

dx, 


denote the deformation gradient, its determinant, and its inverse, respec- 
tively. The material time derivative of the second Piola-Kirchoff stress tensor 
and the Truesdell rate of the Cauchy stress tensor are related in the same 
manner as the stress tensors themselves (see, e.g. [46]) as 


-jyl ( Skl ) — j XK,k Xl,i ( &ki ) • (2.4.10) 

The Truesdell rate equation is derived by transforming (2.4.8) to spatial 
coordinates as 


j X Atk Xbj ( < 5 */ ) = Cabcd x i,c Xj,D , 

and after further manipulation (2.4.2) results with the definition 

1 ^ 

Cijkl = J Xi,A x j,B x k,C x l,D CaBCD ■ 

In the same manner, a constitutive relation suitable for the present 
beam formulation which incorporates the convected frame decomposition of 
the Cauchy stresses is derived. The following relationship similar to (2.4.9) 

dX K dX L 


>KL = J cr 


pq dir dt q 


(2.4.11) 


is defined by transforming the spatial coordinates from the inertial basis 
to the convected basis. The convected frame interpretation of the Truesdell 
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stress rate is derived by taking the material time derivative of the right-hand 
side of (2.4.11) to yield 


D f „ v . dX K dX L ( * o \ 

m {SKL > = } WW ( u) ' 

and the objective convected coordinate stress rate tensor becomes 


(2.4.12) 


°kl = - a il V ki ~ a ki v li + + a ilPik + a ki hi • (2.4.13) 


In the above equation, the notation 

< = T„„^ , h = % :T ki 

denotes the transformed velocity gradient tensor and the angular velocity 
tensor of the convected reference frame respectively. The constitutive law 
is then derived as 


a ij — c ijkt &kl 


(2.4.14) 


where 


" mp 


I (rp du k . rp <^k X 

“ 2 ( Tmk di P + Tpk dU > 


is the symmetric part of the velocity gradient tensor analogous to (2.3.9). 
The constitutive law is rewritten in terms of skew-symmetric tensors as 

&kl "b a ki ( — Pil ) — ( fiki ~ flki ) a il 

= ( Cklmp + Cklmp ) £mp (2.4.15) 


1 (rp d(lk_ _ rp Jfok X 


diii 


Ump 


di 


dir 


where 
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is the skew-symmetric part of the transformed velocity gradient tensor. It 
is noted that the skew symmetric tensor of the above constitutive law rep- 
resents the difference between the convected frame spin tensor and the con- 
vened frame angular velocity vector. It can be shown that this difference 
is approximately the order of magnitude of the incremental rotations of the 
shear S rotation matrix. Under the assumption of small shear strains, this 
skew tensor can be neglected, and the constitutive equation reduces to 

**kl — Cklmp £mp ■ (2.4.16) 

This equation cam then be integrated directly without requiring a transfor- 
mation of basis as 




+ Cklmp A£ mp 


(2.4.17) 


to define the stress update procedure. 

The above expression (2.4.17) is then written in terms of the relevant 
convected frame stress components as 


r i n+1 f*«r _ 

| | — | / + C < Ae^ > (2.4.18) 

l ) ( °ir, J { J 

The constitutive matrix is subsequently given as 


C 


Cnn C1112 

c 1211 C1212 

Cl311 C1312 

Diag ( E + t 7^ 


C^ 1113 
^1213 
Cl313 


G + 2 ff « 


G + ) 


by specializing the definition (2.4.4) to the tensor components representative 
of the beam formulation. Alternative stress-dependent constitutive laws, 
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as for inelastic material behavior, may also be implemented. The stress 
components axe then written as 

cr = < 7 ° + (.2 < 7 t2 + h ° ts . (2.4.19) 

illuminating the dependency on the cross-sectional variables t { . This de- 
pendency is immediately seen from (2.4.18) due to the incremental strains 
which will be of the same form as the virtual strains (2.3.19). Likewise, the 
stress dependent constitutive matrix is decomposed as 

C = C° + li C /j + £3 C /3 (2.4.20) 

where 

C° = C ( a° ) , C * 3 = C ( <T e * ) , C * 3 = C ( a' 3 ) 

are defined from the appropriate stress components of (2.4.19). By combin- 
ing the above with the following incremental analog to the virtual strains 
(2.3.19) 

Ae = A 7 + t? A/c , 

the stress update thus takes the form 


<7° " +1 

= <7° " 

+ 

A< 7 ° 

a' 3 n+1 

= a h n 

+ 

A< 7 ^ 3 

^3 n- hi 

= a * 3 n 

+ 

Aa t3 

Act 0 = C° 

A7 




Act * 3 = C <3 A7 + C° l 2 A/C 

Act ' 3 = C' 3 A7 + C° ij A/c . 


where 
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The evaluation of the strain increments A 7 and A« will be detailed in 
Chapter III. The stress resultant forces and moments per unit length N and 
M defined in (2.3.24) are now given in a more explicit form as 


' Ni ' 




' Mi ' 


a 1 " 1 

h - h ' 

n 2 

► = i 

"h A 

► i 

m 2 

► = < 


d? 12 




t 




-4 * ■ 


2.5 Consistent Linearization 

The derivation of the equations representing the spatial motion of 
a flexible beam is thus concluded. Inherent in these equations is an inertia 
force containing a nonlinear rotary acceleration, and a nonlinear internal 
force expression. To achieve numerical solutions to both static and dynamic 
problems, it is necessary to work with a linearized set of equations repre- 
senting the best linear approximation of the nonlinear equations within a 
small neighborhood of some equilibrium configuration. The term consistent 
linearization refers to a process of achieving this approximation in a manner 
that retains the inherent properties of the formulation. Consistent lineariza- 
tion techniques which account for the effect of finite rotations in deriving 
tangent stiffness matrices have been presented for nonlinear beam formula- 
tions [25-27,31-36,54]. In what follows, an intuitive approach is presented 
to linearize the inertia and internal force terms. 

To this end, the displacements and rotations which inherently define 
the nonlinear functions are represented as a linear perturbation about a 
known equilibrium state. An incremental displacement Ax represents the 
linear portion of any possible perturbation of a known displacement x n and 
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to an alternative displacement x as 

x ~ x n + Ax . (2.5.1) 

The equivalent concept which gives the linear portion of any possible ex- 
cursion between a known orientation represented by the orthogonal rotation 
matrix R" and an alternative rotation R is given as 

R ~ R n + AR (2.5.2) 

= R" + Ad r R" . (2.5.3) 

The definition AR = Aa T R n is analagous to the virtual rotation defini- 
tion <5R = 8a T R of (2.2.9); the components of A a T represent infinitesimal 
angles of rotation about the basis vectors of R n moving the reference frame 
from an unperturbed to a perturbed position [4]. The form of the linearized 
rotation has also been derived using the Frechet derivative concept in [31,34]. 

2.5.1 Tangent Mass 

To obtain the linear portion of the inertia force 


II 

f { 6ut 

Sa T } | 

Mm 1 

(2.5.4) 


= pA il 



(2.5.5) 

Mj 

= J Co 

-j“ UJ J 

U , 

(2.5.6) 


linear perturbations of the translational and rotational acceleration must be 
derived. The former is given as 


u 




il n + A u 


(2.5.7) 
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such that 

M m ~ + pAAu . (2.5.8) 

The latter will be obtained from linear expansions of the angular veloc- 
ity and acceleration in terms of the rotational perturbations about a given 
orientation. For this purpose, (2.5.3) is employed to give 

(5 = -R R r 

~ -j t [ ( I + Aa T ) R" ] R" T ( I + Ad ) . (2.5.9) 

By expanding the above and retaining terms to first order in the rotational 
perturbation, the consistent approximation of the angular velocity tensor 
becomes 


u ~ u> n + Ad + A n Aa - Act £j n . (2.5.10) 

The vector dual to the above skew-symmetric tensor is given as 

u ~ u n + Ad + w"Aa . (2.5.11) 

Linearization of the angular acceleration is derived in the same manner from 
the definition 

<5 = -RR r - RR r , 
and the result is given in vector form as 

ui ~ u n + a + u 1 Aa + u n Aa . (2.5.12) 

From (2.5.11) and (2.5.12), the linear expansion of the rotational accelera- 
tion can then be derived as 

Mj ~ Mj + J A0 + C G A 0 + K c A 0 (2.5.13) 
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where C G and K c represent gyroscopic damping and centrifugal stiffening 
matrices respectively given by 

C G = u n J - 3u n + Ju 11 (2.5.14) 

K c = A n JA n - (Ju n )u n J + Ju 1 ■ (2.5.15) 


2.5.2 Tangent Stiffness 

To obtain the linear portion of the internal force 

5F S = j v T m ,^ dV , (2.5.16) 

linear perturbations of the convected frame transformation matrix, the spa- 
tial derivatives, the stress state, and the deformed volume from a given de- 
formed position must be derived. The perturbation of the rotational tensor 
T{j representing the convected reference frame is given as 


Tij a ( S ik ~ A 0 ik ) T k] n (2.5.17) 


where A/?,* are defined as linearized rotation increments of the convected 
reference frame. To obtain a linear perturbation of the spatial derivative, 
(2.5.1) is employed to give 


dxi 

dx } n 

This is then inverted as 


( + 


dAxi N 

dxf ’ 


(2.5.18) 


dx x n 

dxj 


~ ( S- - gA — ) 

— V °*J a _ n ) 


d Xj n 


(2.5.19) 


as the displacement increments Ax{ are assumed to be small. The spatial 
derivative can thus be approximated as 


dxj n 


dxi dx- n 


dAx t . d 
' Sij dx n ’ 


dxf 


d_ 

dxi 




(2.5.20) 
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The above can then be transformed to convected coordinates as 


a _d 

d(i ~ Tij dXj 

* ( tit - Wit ) 3$ ( S ri - ^ ) 9 


dxf 


dx p n 


(2.5.21) 


by incorporating (2.5.17) and (2.5.20). The final approximation is given as 


4 ~ ( S - k - - V jfi ) ■ (2.5.22) 

by transforming (2.5.21) back to convected coordinates and retaining all 
terms to first order. A consistent perturbation of the deformed volume is 
obtained in a similar manner as 


iV = iet{ £? ) dV " 

= det(Sa + 

(1 + T m ; ) dV« . (2.5.23) 

u Cm 

A consistent perturbation to the Cauchy stress state is defined from the 
constitutive law (2.4.13) and (2.4.14) as 




~ a 




+ A<7,/ 


(2.5.24) 


The proper stress increment is deduced from the time derivative of the con- 
vected Cauchy stress components. This is obtained by rearranging (2.4.14) 
and conceptualizing ( A ) to give 




( C ijkl - SuaJ + 6 u <r jk a + 6 ]t a lk a ) T ks ^ 

+ cr,/A fa + cr u a A(3tj . (2.5.25) 
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The linear perturbation of the internal force can be identified by 
substituting the approximations (2.5.17, 2.5.22 - 2.5.24) into (2.5.16). After 
an extensive set of calculations with some convenient algebraic cancellations, 
the result becomes 


8F S 

~ 8F S + 





■ 
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d Ax. 


Tki 

d8ii a 

dAx s 

- dv . 


/ T ki 

J V 

dil Ck,m r Tma di P 
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di, c ’ , ‘ Tk ' 

dt P 

This is 

further 

simplified as 









8F S ~ 

C JV 

8F S + 

8 K m 

+ 8 K g 


(2.5.26) 
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j ^ kl Cklmp 

temp dv 


(2.5.27) 



8 K g = 

^ 1/ 

f d8ii 

J. % 

a pi 

dAxi 

dv 


(2.5.28) 


in which the material tangent stiffness operator 8 I\ ^ and the the geometric 
tangent stiffness operator SK G have been identified. 


Explicit expressions for the material and geometric tangent stiffness 
operators are given as follows. The material stiffness is simply 

8 K m = J^{8u T 8a T } B t CB j (2.5.29) 

where B is the strain-displacement operator defined in (2.3.28) and 


C = Diag ( EA , GA , GA , GJ , EI2 , EI3 ) 


is the elastic constitutive matrix incorporating the area integration through 
the cross-section variables. To obtain the geometric tangent stiffness, 
(2.5.28) is written in terms of the relevant components as 


r f d8r 

K = l 


d8r dAr / d8r dAr d8r dAr 

d£ di V d£ drj + drj d£ 

( d8r dAr d8r dAr ^ 

v Hz' !k~ + Ik ~dT ' 


/ d8r dAr d8r dAr \ 
v d£ dr] + dr] d£ ) 


(2.5.30) 
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The expressions (2.3.15 - 2.3.17) as well as (2.4.19) are then substituted into 
above and the integration through the cross-sectional area is performed. The 
result becomes 


SK l 


-l 


dlu T „ 3&u 

( ~sr Nl ~sr + 

d6u T T dA/3 d8u T T 
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(2.5.31) 


where the stress dependent matrices introduced in the above are given as 
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Finally, the convected frame virtual rotations Sj3 and their spatial derivatives 
are transformed back to the cross-section reference frame via (2.3.26) or 
(2.3.27). The tangent stiffness matrices K M and I\ G can thus be identified 
from the above derivations such that the lineaxized internal force is given as 


6F S = { Su T Sa T } 5 dt 

~ { 6u T Sa T } | S n + 


[*° + KM i {£}} 

(2.5.32) 
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2.7 Concluding Remarks 

A formulation for the dynamics of flexible beams that admits the 
effects of large rigid-body motions as well as large deformations has been 
presented. The present formulation adopts an inertial reference frame for de- 
scribing the translational motions and a body-fixed frame for the rotational 
motions. As such, the beam inertia operator is identical in form to that of 
rigid body dynamics and the formulation is easily adapted into the method- 
ologies of general multibody dynamics. To effect the kinematic description, 
nonlinear strain theories are incorporated from the outset which account for 
both finite rotations and large deformations. The present nonlinear strain 
measures are deduced from the form of the virtual work of the internal force 
of the beam as parameterized by the kinematic definition. An rate- type con- 
stitutive law is employed to complete the formulation. A key feature of the 
formulation is the use of a convected reference for the stress representation. 
This leads to conceptual simplifications which enable a Cauchy stress repre- 
sentation of the beam deformation as opposed to the classical Piola-Kirchoff 
stress representations typically used in finite-deformation analysis. A con- 
sistent linearization of the equations are presented for reference in future 
chapters in which the tangent stiffness matrices are inherently symmetric as 
a consequence of the Cauchy stress representation. 

The next chapter discusses the discretization of the equations of 
motion and the numerical procedures developed for the computation of the 
discrete internal force expression. 



CHAPTER III 


COMPUTATIONAL TREATMENT OF THE 
DISCRETE INTERNAL FORCE 

3.1 Introduction 

The previous chapter has presented a formalism for modeling the 
spatial dynamics of flexible beams. The formalism accounts for both the 
finite rotations and the finite deformations of a beam component. The dy- 
namics of the beam motion are described using an inertial reference for 
translational displacements and a body-fixed reference for rotational quan- 
tities. To effect this description, finite strain rod theories are defined. 

The numerical treatment of the formalism is now considered in detail 
in the next two chapters. The present chapter is devoted to the numerical 
treatment of the internal force operator, and the next chapter discusses the 
dynamic solution procedures. For the effective treatment of the internal 
force, a computational invariance to rigid body motion must be satisfied 
within a discrete beam model to achieve realistic static and dynamic simu- 
lations. As the degrees of freedom contain information of both rigid motion 
and strain deformation, the internal force computation depends on a ju- 
dicious procedure which filters out the rigid motions implicitly embedded 
within the variables. This internal force computation can then be interfaced 
with multibody dynamics solution procedures, to be discussed in Chapter 
IV, with confidence that spurious strains are not computationally generated 
in the process. For this purpose, a unique procedure for the computation 
of the internal force has been developed by exploiting characteristics of the 
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strain formulation as well as finite rotation theory. 

The finite element method is employed to give a finite-dimensional 
representation of the continuum beam model. The variational form of the 
partial differential equations presented in Section 2.3 provides the basis for 
a displacement-based finite element methodology [64], resulting in a set of 
ordinary differential equations to be integrated by procedures discussed in 
Chapter IV. As such, the beam is represented as an assemblage of elements 
interconnected at nodal points on the element boundaries, and the degrees 
of freedom of the beam become the position and orientation coordinates of 
plane sections at the nodal points. 

This procedure is a unique contribution to the computational re- 
search in the area of flexible multibody dynamics. Previous inertial refer- 
ence beam formulations are based on total strain measures which are defined 
with respect to constant material coordinates of an initial configuration [25- 
27,31-36]. The invariance of the particular strain computation has not been 
specifically proven within these works. The present work is based on an 
incremental formulation, thus providing a more natural interpretation of 
large deformations by retaining the history of the motion in determining 
the current stress state. This formulation differs from the previous works 
in which the stress state is determined by comparing the current shape of 
the continuum to its original reference configuration. In a manner similar 
to the previous works, the present formulation employs non-commutative 
orthogonal matrices for the general representation of finite rotations. A 
complete freedom is thus allowed in choosing the parametrization of the 
orthogonal transformation which achieves the most effective computational 
scheme. The present work exploits this freedom by incorporating various 
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parameterizations to derive accurate methods for the computation of strain 
increments within the stress update procedure. 

The stress within a beam finite element is computed from the nodal 
translational and rotational quantities as follows. Previous inertial reference 
beam formulations are based on total strain measures which are defined with 
respect to constant material coordinates of an initial configuration [25-27,31- 
36], As such, the stress state is determined by comparing the current shape 
of the continuum to its original reference configuration in these works. The 
present work is based on an incremental formulation, thus providing a more 
natural interpretation of large deformations by retaining the history of the 
motion in determining the current stress state. As discussed in Section 2.4, 
the present formulation adopts an incremental procedure based on rate- 
type constitutive laws in which the current stress state is updated from the 
past stress state via a strain increment and appropriate constitutive ma- 
trix. A computational procedure has been designed from this incremental 
interpretation of the continuum-based formulation such that the computed 
finite strain increments are invariant to arbitrary rigid body motions implic- 
itly contained within the finite displacements. To achieve this invariance, 
various rotational parameterizations axe incorporated in deriving accurate 
methods for the computation of strain increments within the stress update 
procedure. As will be shown, the success of the computation hinges on 
a proper extraction of the rotation increments from the various rotational 
matrices embedded in the formulation. Given the strain increments, the 
elemental internal force is then formed via the multiplication of a discrete 
strain-displacement operator with the updated stress state. 

The computation of the elemental internal force thus involves the 
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proper treatment of large rotations. A great deal of literature exists on 
the subject of finite rotations. An extensive derivation of the kinematics 
of large rotations is presented in [55], and an overview of different methods 
for parameterizing an arbitrary rotational orientation is given in [56]. A 
brief overview of the various parameterizations employed in previous works 
involving large rotations is given as follows. Euler angles or Bryant angles, 
which are based on successive angular rotations about a predetermined set 
of axes, were historically used for specialized rigid body systems and at- 
titude dynamics [47-48]. To alleviate problems associated with the non- 
commutative nature of this type of definition, commutative semitangential 
rotations were derived and applied to a large displacement-small strain struc- 
tural theory [54]. The rotational vector, which is based on the Euler-Chasles 
description of a single rotation about an appropriate axis [55], has been used 
within fully nonlinear beam and shell theories [31-34,57]. While the angle- 
based geometrical interpretation of these parameterizations is recognized, 
the rotational transformation matrix and subsequent angular velocity and 
acceleration vectors are complicated trigonometric functions of these rota- 
tion parameters. To avoid the complications of these nonlinear functions, 
alternative parameterizations were developed such that the rotational ma- 
trix, the angular velocity vector, and the angular acceleration vector become 
algebraic functions of the rotation parameters. Examples of these alternative 
parameterizations include the Rodrigues parameters [28,55], the conformal 
rotation vector [36,58], and the Euler parameters [59-63]. In general, the 
algebraic nature of these parameterizations is exploited when developing 
computational strategies for general finite element and articulated system 
models. The present formulation incorporates both the Rodrigues parame- 
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ter and the Euler parameter representations in developing an efficient com- 
putation of the internal force which remains invariant to arbitrary rigid 
body motions. The Euler parameter representation is further exploited in 
Chapter IV for the integration of general spatial kinematic systems. 

The computational procedure developed in this chapter is compared 
to previous works through an analysis of static problems. The nonlinear 
static equilibrium equations are solved in an iterative manner using the 
Newton-Rhapson technique. The procedure is included for verification pur- 
poses as the present formulation is intended more for use with multibody 
dynamic applications rather than to compete with existing large rotation 
and deformation beam finite elements used within static analysis. 

The rest of the chapter is organized as follows. Section 3.2 gives 
an overview of the rotational vector and Euler parameter representations 
of finite rotations. The constant-strain finite element discretization of the 
equations of motion is presented in Section 3.3. The algorithmic treatment 
of the internal force is detailed in Section 3.4, and the static solution verifi- 
cation is given in Section 3.5. 

3.2 Finite Rotation Representations 

The Euler- Chasles theorem states that any finite rotation can be 
uniquely represented with a rotation angle 6 and a rotation axis n [47-48]. 
This representation is shown below in Figure 3.1. 

From this picture, we see a new vector r r is obtained by rotating a given 
vector r by an angle 6 about an axis n. The mathematical description of 
this process is given by [47-48] 


r = 


R r 


(3.2.1) 
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Figure 3.1 Euler-Chasles rotation representation 


where R denotes the dydadic form of the rotational operator given as 

R ( n, 6 ) = nn + cos# (I — nn ) + sin#(Ixn) (3.2.2) 

with I representing the (3 x 3) identity operator. The above definition is 
consistent with that of the angular velocity given in (2.2.9). The matrix 
representation of the above definition is given as 


R = / + + L e e 


6 


# 2 


(3.2.3) 


In this representation, the parameter 

0 = n# 

is termed the rotational vector, and 0 is the skew-symmetric dual matrix 

0 = 


0 1 1 
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02 0 —0i 

-02 01 0 
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defined from the three components of the rotational vector. The matrix 
(3.2.3) is equivalent to the following series expansion [55] 


R 


/ + 0 T + 





+ ••• + 



+ ••• 


The above representation is precisely the matrix exponential 


R = exp ( 0 T ) , (3.2.4) 

and thus the rotation matrix is a nonlinear exponential function of the 
rotational vector components. 

An incremental rotation which relates an existing body-frame basis 
b n to a new basis b n_ *’ 1 is defined as follows [56]. Given that the existing 
basis is defined with respect to an inertial basis e by the rotation matrix 
R n as b" = R n e, and likewise a new basis is defined as b n+1 = R n+1 e, 
the rotation matrix R”"*" 1 is obtained from the product of an incremental 
rotation matrix R and the existing rotation matrix R n . This matrix product 
can be performed as 

R n+1 = R(/> R" (3.2.5) 

in which the incremental rotation operator R(i), termed a spatial rotation, 
is applied to the existing body frame b". A material rotation can also 
be defined in which an incremental operator R( r ) is applied to the inertial 
reference frame as 

R n+1 = R n R (r) - (3.2.6) 

The rotational vectors 0 and $ that correspond to the spatial or material 
incremental rotations R(/) and R(r)? respectively, can be introduced as 

R(0 = e® T R (r) = e* T 


(3.2.7) 
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The two rotational vectors are related via 0 = R n $ . 

While the definition of the rotational vector is straightforward, the 
drawback to this parameterization is the nonlinear expression for the rota- 
tional matrix. To avoid the nonlinearity seen in (3.2.3) or (3.2.4), alternative 
parameterizations have been developed which lead to algebraic expressions 
of the rotational matrix. For example, the Rodrigues parameters are defined 
as 

Q 

0< = 0 t n , 8 t = tan - (3.2.8) 

such that the rotational matrix is given as a function of these parameters as 

R = 1 + 1 l 0 2 (©f + ©*©<) • (3.2.9) 

The drawback to the above parameterization is that it becomes singular 
whenever the magnitude of the rotation becomes an odd multiple of 7 r. In 
general, a singularity condition is encountered within any three-parameter 
representation of orthonormal rotation matrices [65]. An alternative non- 
minimal parameterization, termed the Euler parameters, remains singularity 
free throughout all possible orientations. The Euler parameters are defined 
as 

0 0 

q 0 = cos - , q = sin - n , (3.2.10) 

and the four parameters are subject to the constraint 

9o + q T q = 1 • (3.2.11) 

The rotation matrix is given as a function of the Euler parameters as 

R = I + 2? 0 q T + 2qq . (3.2.12) 
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Given these examples, for computational purposes it is necessary 
to extract the various rotational parameterizations from given rotational 
matrices. If an inverse exists, the rotational vector can be obtained from 
(3.2.4) as 

Q T = log e R = -sin- 1 I R° (3.2.13) 

T Z 

where 

R a = R - R t (3.2.14) 


is formed from the anti-symmetric part of R and t is the norm of the dual 
vector to R° [66]. For the Rodrigues parameters, the inverse of (3.2.8) is 
the following algebraic expression [28] 


Qj = 


R“ 


1 + tr R 


(3.2.15) 


For the Euler parameters, a very robust method involving only algebraic 
computations has been derived in [67] to extract the Euler parameters from 
the rotation matrix. 

To complete this introduction to finite rotations and their param- 
eterizations, it remains to express the angular variation, angular velocity, 
and curvature vectors as functions of the given rotational parameterizations. 
The body-fixed components of the skew-symmetric tensors representing the 
angular variation, angular velocity, and curvature are defined as 

Sa T = 6R R r . <i T = R R r , k T = ® R T , 


respectively. These definitions constitute spatial representations; the corre- 
sponding material representation is given as 



R t Sa T R 


R r 6R 
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for the angular variation with corresponding definitions for the material 
angular velocity and curvature tensors. From the above definitions, it is 
seen that the derivative of the rotation matrix is necessary to obtain explicit 
expressions for these functions in terms of a desired parameterization. 

For the rotational vector parameterization, the derivative of the 
exponential map is derived in [31,66]. The spatial angular variation is then 
given as a function of the spatial rotational vector and its derivative as 

~ ^3 (© ' ^ 0 ) 0 + ~ Q — ^0 + 02 ( ~ ^00 ) • 

(3.2.16) 

The equivalent vector expression for the above skew-symmetric tensor is 
given as 

6a = T r (0 ) 60 (3.2.17) 

where 


Tr(0) = 


9 — sind 

¥ 


0 0 r + 


sind 


I + 


1 — cosQ ~ T 
02 0 


The material components of the angular variation are related to the material- 
based rotational vector definition in the same manner as 6a e = T r ($ ) 0$ 
[56]. Analogous expressions exist for the angular velocity and curvature ten- 
sors. 

For the Rodrigues parameters, the functional relation 


6a = T, (0 t ) 00, 


(3.2.18) 


results [56], where 


T, (0, ) = 


1 + 6 } 


(i + $n 


(3.2.19) 
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The identical relation holds for the material components of the angular 
variation using the material components of the Rodrigues parameters. 

The corresponding functioned relation in terms of the Euler param- 
eters is given as 



?o qT 1 f S( lo \ 

-q ?oi - q J l £q / 


(3.2.20) 


This expression, can be inverted to give the Euler parameter derivatives in 
terms of the body-frame angular variation components as 





in which the constraint conditions 


( 3 . 2 . 21 ) 


+ q r q = 1 , <$<Zo9o d" ^ qT q — P (3.2.22) 

are embedded. The corresponding inversions of (3.2.17) and (3.2.18) have 
been derived in [66]. 

The Euler parameters do not posses any singularity limitation 
within their definition, nor within the forward or inverse transformations 
between the angular and parameter variation. As such, the Euler param- 
eters are the only set that allows the treatment of rotations of arbitrary 
magnitude without resorting to any special precautions. In contrast, the 
Rodrigues parameters possess an unavoidable singularity in their definition 
[28], and the rotational vector parameterization possess a “hole in differen- 
tiability” in which the transformation matrix of (3.2.17) becomes singular 
[34]. Thus, the Euler parameter representation of the angular velocity and 
curvature has been chosen in the present work in developing the computa- 
tional procedures to be discussed shortly. For the representation of incre- 
mental rotations, both the rotational vector and the Rodrigues parameters 
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are introduced where appropriate, again with the motivation of achieving 
effective computations dictating the particular choice of the parameteriza- 
tion. These concepts are referred to and further illuminated in the following 
sections. 

3.3 Finite Element Discretization 

A finite element representation of the continuous beam model is 
obtained from the variational form of the partial differential equations pre- 
sented in Section 2.3. The present study is based on the use of linear shape 
functions in approximating the displacement field within a beam finite ele- 
ment [64], viz., 

npe 

« = £>, ( X ) Ui ( t ) (3.3.1) 

7=1 

where Nj denotes the spatial linear shape functions, u j represents the de- 
grees of freedom at the element nodes, and npe denotes the number of nodes 
per element. In the present beam formulation, the rotational variables are 
treated independently from the translational variables. As such, only Co 
displacement continuity across element boundaries is required for the shape 
function approximations. Thus, the linear interpolation functions result in 
a consistent approximation which leads to a constant strain finite element. 

As defined above in (3.3.1), the interpolation of the translational 
quantities is straightforward as the components are expressed with respect 
to a common reference frame. However, the interpolation of the rotational 
variations requires special consideration. The virtual rotation components 
8a used within the present formulation are defined with respect to a body- 
fixed reference frame which is continuously changing along the beam neutral 
axis. As such, these virtual rotation components cannot be directly inter- 
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polated, but rather must be transformed to a parameterization defined with 
respect to a fixed configuration prior to the interpolation. To this end, the 
virtual rotations can be parameterized in terms of the rotational vector via 
(3.2.17), where 0 is the rotational vector orienting the moving basis from a 
common reference as the inertial frame. The components of the rotational 
vector are defined with respect to a fixed coordinate system as 

0 = 8 n = 9\ e, , 

and hence the variations of the rotational vector can be interpolated as 

npe 

<50 = ^ Ni 60/ . 

i=i 

From the above, the virtual rotations 6a can then be interpolated as 
<5a = T(0) <50 


npe 

~ T(0) Ni 60/ 

1=1 

npe 

= Ni T(0) T _1 (0/) 6a/ . 

i=\ 

If 0 =» 0/, the product T(0) T _1 (0/) =» I [34]. Thus, if the rela- 
tive change of the cross-section orientation along the finite element remains 
small, the virtual rotations are approximated in the same manner as the 
displacements as 

npe 

8a ~ ^ iVj 8ai . (3.3.2) 

1=1 

The angular velocity and angular acceleration are then also interpolated as 

npe npe 

ui ~ N/ u>/ , cj ~ N/ u>/ 

1=1 i=i 


(3.3.3) 
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From this introduction to the finite element approximations, a discrete form 
of the present beam formulation is readily obtained as follows. 

3.3.1 Discrete Equations of Motion 

A discrete form of the variational equations (2.3.1) is obtained from 
the finite element approximations 

npe. npe 

8u = J2 N I 6u i » Sa = Y, Ni t<*i (3.3.4) 

i=i i=i 

for the virtual displacements and virtual rotations, respectively. As standard 
linear shape functions are chosen for Ni, uj and 8a j represent the quantities 
at the two element nodes. From the above interpolations, a discrete form of 
the inertia, interned, and external force operators presented in Section 2.3 
axe obtained for a single finite element. Finite element assembly procedures 
are then applied to the element operators to achieve discrete variational 
equations in terms of the nodal degrees of freedom of the beam configuration. 
To this end, the discrete form of the elemental operators are given as follows. 

The internal force operator of an individual finite element is ob- 
tained in terms of the virtual displacements and rotations at the element 
nodes from (2.3.25) and (3.3.4) as 


"PC 

6F S = Y, {<*«/ 6a i) 

i=i •*£ 


(3.3.5) 

[s £ ,] = 

rr dN, 

A 

*i Nj S J 

(3.3.6) 


0 

Sj(ksNj + f ). 


By introducing linear shape functions into the above and integrating across 

the element length, the above equation 

can be written as 


8F S = {<5ui 8u2 <5<*i 


(3.3.7) 
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where 

[b e ] = 

' -*T IT |fi Sf i*J sr 

0 0 Sf ( J*s - }I ) S* ( l*s + ll ) . , 

(3.3.8) 

In deriving the above expression, an underintegration of the rotational vari- 
ables associated with transverse shear strain has been performed to avoid 
the locking phenomena associated with the Timoshenko beam formulations. 
The convected frame T matrix, body frame curvature tensor k s , element 
neutral-axis length £, and resultant stresses 1V 7 and M K are constant quan- 
tities over the element domain, while the relative cross-section deformation 
S matrices axe nodal quantities as denoted by their subscripts. 

In the same manner, the element external force operator is obtained 
from (2.3.29) and (3.3.4) as 

6F b = £ {<«, to,} , ff = J , 

7=1 (3.3.9) 

and the traction operator is implemented as boundary conditions on the 
nodes. Likewise, the element inertia operator is obtained from the (2.3.4), 
(3.3.3), and (3.3.4) as 

npe npe *2 

6F- = { 6uJ f>AM B , K ig- } + 

7=1 K=\ 

npe npe , 6 . . 

£E { ScJ.J, } + 

7=1 K=l 
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where 


npe 

S>Fo E / 


/=1 


( 3 . 3 . 10 ) 


M 


E 

IK 


t, 


N, N k d( 


d e M, 



represent the element mass matrix and nonlinear angular acceleration vec- 
tor. The former will be evaluated as a standard lumped mass matrix for the 
computational efficiency of explicit integration techniques to be described 
in Chapter IV, and the latter will be evaluated in a consistent manner by 
averaging the element nodal angular velocities. 

When standard finite element assembly procedures axe applied to 
the element operators presented above, one obtains the discrete variational 
equations 


{ 6u/ Sa d T } 


- 

r m d u d ' 


0 

-1 

< 

1 1 

> + < 

( + \ 

► 

■ 

[ Jd J'd J 


{ D d (u) ) l S,$ J 

- 


= {Su d T Sa d T } |^| ( 3 . 3 . 11 ) 

where 6u d and 6a d now represent the nodal virtual displacements and rota- 
tions of the assembled finite elements modeling the beam configuration. The 
following notation has also been introduced in the above equation; m d and 
Jd represent the assembled mass and inertia matrices; iid and u d represent 
the nodal acceleration vectors; D d { u>) represents the assembled nonlinear ac- 
celeration, and S/' b and f d ' b represent the assembled internal and external 
force vectors partitioned into translational and rotational parts, respectively. 
As Su d and 6a d represent arbitrary independent variations, the final discrete 
equations of motion are given as 
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rrid 0 

0 J d 


u d 


+ 


0 

DdH 


+ < 


< S3 1 

l SJ J 


> = 


ft J 


(3.3.12) 

These equations can be specialized to the case of static equilibrium as 


S d — fd • (3-3.13) 

The static equilibrium equations provide the basis for the remainder of the 
chapter which is devoted to the computation of the internal force. The 
extension of the computational procedure to include the dynamics of the 
inertia operator is presented in Chapter IV. 

To present the algorithmic treatment of the discrete internal force 
S d , the rest of the chapter is organized as follows. The evaluation of (3.3.S) 
is detailed by presenting the computation of the discrete rotation matrices in 
Section 3.3.2 and the computation of the element curvature in Section 3.3.3. 
The computation of the resultant element stresses iV T and M K , necessary 
to complete (3.3.7), is then presented in Section 3.4. Finally, the numerical 
solution of the static equations (3.3.13) is presented in Section 3.5. 


3.3.2 Discrete Rotation Matrix Computations 

The rotational matrices representing the various reference frames 
introduced in the present formulation are computed as follows. The compu- 
tational procedure for the solution of the discrete equations of motion will 
be discussed in Chapter IV. From this solution procedure, the orientation 
of the cross-section and the displacement of the neutral-axis at the finite 
element nodes are output as known quantities. The matrices R,, repre- 
senting the cross-sectional orientation at each element node, are computed 
directly from the rotational parameterization of the given solution method. 
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The rigid motion of -the convected reference frame a of a given finite ele- 
ment is then deduced as follows. The translational displacements axe used 
to compute the tangent to the neutral axis ax ; for the linear finite element 
interpolations the tangent is simply the straight fine connecting the two el- 
ement nodes. The remaining two basis a 2 and a 3 of the convected reference 
frame axe then constructed from the tangent aj and the cross-sectional ori- 
entation Rx . The a 2 vector is defined as the cross product of a x with the 
b 3 axis of Rx, and the a 3 axis is then defined to complete the right-hand 
coordinate system. The computed axes { ax,a 2 ,a 3 }, as shown below in 
Figure 3.2, define the rows of the T matrix. 



Figure 3.2 Convected Reference Frame 


The rotation matrices S, , defined at each element node as the relative 
difference between the element convected frame and the nodal body frames, 
axe thus 


S, 


R, T r , 


i = 1,2 . 


(3.3.14) 
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The computation is an approximation applicable for moderate strains such 
that the Sj matrices contain information solely due to transverse shear and 
torsional deformations [28]. 


3.3.3 Discrete Curvature Computation 

To complete the evaluation of (3.3.8), the discrete computation of 
the element curvature is presented as follows. The body frame components 
of the curvature tensor k^ can be equivalently defined using the R or S 
matrices as 


0 «3 — «2 

— «3 0 «1 

«2 — «1 0 



d e R 

~W 


K 1 


(3.3.15) 


since by construction the convected frame T matrix is constant along the 
element domain where the differentiation is performed. Thus the curvature, 
which was originally defined from the shear matrix in (2.3.27), can be di- 
rectly computed from the nodal rotational variables of R denoting the total 
cross-section orientation. This is advantageous as the Euler parameters of 
R axe directly available from the generalized coordinate integrator discussed 
in Chapter IV, whereas additional computations axe required to obtain the 
Euler parameters corresponding to the S matrix. 

The Euler parameter representation of finite rotations can be ex- 
ploited to achieve a simple and robust computation of the element curva- 
ture. The Euler parameter - curvature relation, introduced in (3.2.20), is 
given as 


= | 


K = 



(3.3.16) 
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where 


E ( q ) = 2 f 50 qT 1 = I ft \ 

1 ' -q *>i-q J d( l If / 

The above is evaluated from the Euler parameters of the element nodes by 
using the normalized average 


_ I ( 9 * + & ) 
q ‘ ~ II 5 ( 1. + ) II 

for the matrix E and the approximation 


dq 

di 


J ( 92 - qi 


) 


(3.3.17) 


(3.3.18) 


for the derivative. Within the geometric context of rotations, the Euler 
parameters q a of (3.3.17) correspond to an orientation averaged from the two 
nodes. The computation of (3.3.16) from (3.3.17) and (3.3.18) is consistent 
in that the constraint conditions 


9o 2 + q T q = 1 , ^-?0 + q = 0 

are identically satisfied in the discrete sense. Due to the singularity free 
nature of Euler parameters, the computation is valid for any size nodal 
rotation. The use of an alternative rotational parameterization for the com- 
putation of the element curvature requires special treatment to avoid the 
inherent singularities within the representation [34]. 

To show the robustness of the averaging scheme adopted herein, 
we offer the following analysis. By examining the rigid rotation of a finite 
element in a state of constant curvature, it can be shown the computa- 
tional strategy effectively filters out the large rigid-body motions implicitly 
contained within the nodal Euler parameters. For this purpose, the total 
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orientation of the nodal cross-sections are described by separating from the 
outset the effects due to rigid motion and that due to deformation. The 
orientation of the convected element frame is characterized by a rotation 
of an angle 4> about an axis n„ from the inertial reference frame, and the 
relative nodal cross-section orientations are characterized by a rotation from 
the convected frame of angles — r and r about axis n& for nodes 1 and 2, 
respectively, as shown below in Figure 3.3. 



Figure 3.3 Pure Bending of Beam Element 


From these definitions, the Euler parameters q Tl and q r2 representing the 
relative nodal orientations and q a representing the convected orientation are 
given as 

f cosf 1 
\ sin j n b J 

The Euler parameters q\ and <72 designating the total cross section orien- 
tation of the two nodes due to these combined effects can be obtained by 
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applying the quaternion product rule [68] to the above definitions to yield 

_ f cos | cos y -f n„ • n b sin f sin y 1 

\ — cos ^ sin y n t + cos y sin £ n„ — sin y sin | n a x nj / 


cos ^ cos y — n a • nj sin ^ sin y 1 

cos | sin y rift + cos y sin ^ n„ + sin y sin | n a x nj / 


The average of the above nodal parameters is simply 

i(». + ») = { co : fco ^ } 

2 [ cos j sin f n a J 

the norm of which is cos It is seen the normalized average precisely rep- 
resents the orientation of the cross-section at mid-element which contains 


no relative deformation as equal but opposite relative rotations were applied 
to the ends of the element. The normalized average of the nodal Euler pa- 
rameters thus corresponds to an orientation geometrically midway between 
that of the two nodes. It is then easily shown from the constructed parame- 
ters that the computed curvature for the pure bending of the finite element 
becomes 

4 . r 

* = l sm 2 nb 

as opposed to the true curvature of f rnj. It is noted that the above compu- 
tation retains only the rotation parameters r originally defined relative to 
the rigid body orientation, and is thus invariant to any rigid body motions. 
The approximation is valid for small relative rotations of the cross-sections. 
For instances when the validity of the approximation is challenged, an incre- 
mental curvature computation can be made as discussed in the next section, 
from which the total curvature is obtained from an appropriate update pro- 


cedure. 
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3.4 Incremental Strain Computations 

To complete the algorithmic treatment of the nonlinear stiffness op- 
erator, it remains to discuss the computation of the element stresses. The 
stress update procedure, discussed in Section 2.4, requires an effective com- 
putation of the increment of strain between a past converged solution and a 
current solution. A specific computational strategy has been developed for 
use with this incremental interpretation of the continuum-based formulation. 
The fund amen tal property motivating this development is the condition that 
the finite strain increments be computed from the nodal degrees of freedom 
in a mann er that is invariant to any arbitrary rigid body motions contained 
within these variables. Given an accurate computation of the incremental 
strains, the stress update procedure directly follows as discussed in Section 
2.4. 

The incremental strains are defined to approximate the time inte- 
gration of the strain rate tensor over a given time interval. To effect this 
approximation within the present virtual work formulation, the kinemati- 
cally admissible virtual displacements and rotations within the virtual strain 
displacement relations axe defined to coincide with the actual displacements 
and rotations occuring within a finite increment of time. The strain incre- 
ments are thus defined in the same manner as (2.3.20) and (2.3.21) as 


A7 


A K 


dAu 

~w 


+ 


d A/3 


0 

-Aft 

A& 


where Au and A/? are finite displacement and rotation increments, respec- 
tively. For computational purposes, it becomes necessary to separate the 
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rotation due to rigid body motion from that due to deformation within the 
convected frame components of the incremental rotations A/3. The con- 
vected frame virtual rotations 6/3 can be decomposed as 


6$ = 6r a + 6<p (3.4.1) 

where 6(p and 6r a are the convected frame components of the rigid body 
and deformation virtual rotations, respectively. These virtual rotations are 
defined as 

hq> T = 6TT t , 6fJ = S T 6S . (3.4.2) 

The definitions (3.4.1) and (3.4.2) have been derived from the identity 


R = S T , tfR = £S T + S ST 
and the following definitions 


6a T = 6RR r , 6f3 T = S T Sa T S 
6f T = 6S S T , SfJ = S T 8r T S . 


With the decomposition (3.4.1) interpreted incrementally, the membrane 
and transverse shear strain increments are then given by 


x dAu 


f ? \ 

i 9 ) 


^ = T ee 

+ < 

1 —A<p3 > + 

l &P 2 ) 

£ £ 
<!<<] 

(3.4.3) 


Likewise, the incremental curvature representing the torsion and bending 
strains is given by 


as the incremental rotations 
over the element length. 


A k = 


d A r a 

d( 


(3.4.4) 


A <p defined from the T matrix are constant 
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A proper definition and subsequent computation of the finite dis- 
placement and finite rotation increments introduced in these incremental 
strain definitions is essential. The incremental translations are defined as 
the difference between a current and past displacement configuration as 

An = u n+1 — u” . (3.4.5) 

The incremental rotations are interpreted from properly defined rotation 
matrices which represent the difference between current and past orienta- 
tions. It is important to note that the angular variation definitions (3.4.2) 
are consistent with the following spatial and material rotation updates 

T n+1 = AT T n , S n+1 = S n AS , (3.4.6) 

respectively. In the above relations, the matrix AT represents the relative 
orientation between the current and past convected reference frames, and 
AS represents the relative orientation between the current deformation ma- 
trix S n+1 and the past deformation matrix S n rigidly rotated to the current 
convected reference frame. The incremental rotations of the rigid body mo- 
tion, A <p, are then defined as the rotational vector representation of the 
rotation matrix AT, while the incremental rotations of the deformational 
motion, A r a , are defined as the Rodrigues parameter representation of the 
matrix AS. The latter choice is possible as it is highly improbable that a 
singularity is encountered due to large rotational deformations. The differ- 
ent parameterizations are chosen such that objective computations of the 
incremental strains (3.4.3) and (3.4.4) are achieved. 
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3.4.1 Membrane Strain Increments 

To achieve an objective computation, the first two terms of (3.4.3), 

given as 


A rr. ^Au 

A71 = T W + 



(3.4.7) 


must be computed such that rotation increment A <p compensates for the 
rigid motion contained within the displacement increment A u. This is ac- 
complished by defining the skew-symmetric matrix containing the convected 
frame rotation increments as 


A$ t = ( T" +1 - T" ) T n+ 2 T . (3.4.8) 

This definition is a varied form of the approximation 

T n+1 ~ ( / + A$ T ) T n 

for the rotational vector Ay? which parameterizes AT as 

T n+1 = e^ T T” . (3.4.9) 

The reference frame T n+ 2 introduced in (3.4.8) is defined to be geometri- 
cally midway between the current and past convected reference frames as 

T n+ 5 = AT n+ * T" (3.4.10) 

AT" + i = exp ( ) . 

This mid-configuration matrix is necessary in (3.4.8) such that a skew sym- 
metric matrix results for A <y3 r ; this follows from (3.4.8) and (3.4.10) as 

A(f T = T n+1 T n+ ? r — T n T n+ ^ T 
= AT n+ ! - AT n+ ? T 
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To incorporate the computation (3.4.8) within (3.4.7), the matrix T n+5 
must be used within the first term of (3.4.7) to rotate the translational 
derivatives. The reasoning for this is given as follows as well as examples 
exhibiting the invariance of the computation. 

To illustrate that the preceding computation preserves rigid mo- 


tions, let us consider a beam element which is rigidly rotated such that 
A 71 = 0. The nodal incremental translations for the rigid rotation of a 

beam element pinned at the first node becomes 
Aui = 0 



Au 2 = l ( ^" +1 - V 1 ) , u = 

where t ^ corresponds to the first row of the T matrix containing the direc- 
tion cosines of the rotation and t is the length of the beam element. The 


translational derivative within (3.4.7) is then evaluated as 

dAu 


dt 


_ . n+1 x n 

- U ~ *€ 


(3.4.11) 


and the rotation increments are obtained from (3.4.8) as 


| -A<p 3 \ = ~ T n+ * ( tf +1 ~ U U ) ■ (3.4.12) 

l Av?2 J 

With this information, equation (3.4.7) is evaluated as 

A 7l = T-+5 ( ^ - t f " +1 + < f " ) (3.4.13) 

by combining ( 3 . 4 . 12 ) with the translational derivative term rotated by 
T n+ 5. It is then easily seen from (3.4.11) that this particular computa- 
tion leads to 


A71 = T n+ * ( ^ n+1 - if - ^ n+1 + ) 


0 
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as desired. 

As in the preceding example exhibiting the rigid rotation of a beam 
element, an accurate computation of the neutral-axis stretch which is invari- 
ant to arbitrary rigid rotations can also be accomplished. For this case, the 
incremental displacements for an arbitrary rotation and stretch are given by 

Aui = 0 

Au 2 = ((/ + <*) * f " +1 - etf) 

where d represents a stretch relative to the element length i at a past con- 
figuration. The rotational expression (3.4.12) remains valid, and the trans- 
lational derivative is now evaluated as 

{ (* + d) V +1 - it,’-} 

where £* corresponds to the deformed element length at a given configura- 
tion. The above is substituted into (3.4.13) to give 

*71 = T ” +i ( e -jr- <«" +1 - J. V - V + ‘ + V ) (3.4.14) 

for this particular example. From the above, it is immediately seen 
that if the length derivative is evaluated at the current configuration as 
£* = £ t + d , then 


A71 = 


# T ' n +2 + n 

£ + d 1 *€ 


By incorporating (3.4.10), this is equivalent to 


A71 = -7 — — : AT n+ i T n t f n 

£ + d * 


£ + d 


AT n+ i } 


r 1 
0 
0 
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To achieve a computation representative of a pure stretch in which the 
two transverse shear components are identically equal to zero, the above 

, ±T 

computation can be premultiplied by AT” 1 to give 



Thus the general strain computation, obtained by premultiplying (3.4.13) 
by AT n+ s and using the relation 


AT 


/jin+1 



becomes 

r 1 - AT n 1 

A7, = T" ^ + ( AT,, j . (3.4.15) 

It is noted that this computation does not require an explicit evaluation of 
T n+ 2. In effect, the convected frame strain components were evaluated at 
the configuration T n+ s as at this location the computation is invariant to 
rigid motion. The mid-configuration computation is then transformed to 
the current configuration T n+1 to achieve the proper stretch evaluation. 

The same type of result can be achieved in a slightly varied manner 
as follows. The length derivative £* in (3.4.14) can also be evaluated at the 
past configuration as £* = i. Equation (3.4.14) then yields 

A 7 , = . 

This expression must be premultiplied by AT"'*’ 5 to yield a computation 
representative of a pure stretch as 
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Premultiplication of (3.4.13) by AT n+ ^ gives the alternative generic com- 
putation 


A71 = T n+1 


dAu 

&F” 


+ < 



(3.4.16) 


Either evaluation, (3.4.15) or (3.4.16), yields a membrane strain 
computation which is invariant to arbitrary rigid body motions. The con- 
vected frame strain components were evaluated at configuration T n+ 5 in 
(3.4.13) to filter out the rigid motion; this mid-configuration computation 
was then transformed to an appropriate configuration, T n or T n+1 , such 
that the membrane strain component is non-zero and the transverse shear 
strain components are identically equal to zero. To effect a strain evalu- 
ation at the mid- configuration as dictated by the stress-update procedure, 
the computations (3.4.15) or (3.4.16) need only be multiplied by the factor 

i* 

\ ( l n + /"+ 1 ) . 


3.4.2 Transverse Shear and Curvature Strain Increments 

Given the computation of the incremental membrane strain involv- 
ing the displacement and rigid rotation increments, it remains to compute 
the incremental transverse shear and curvature strain increments involving 
deformational rotation increments. To this end, the strain terms 


A72 = 


-At, 

At, 


a 3 


a 2 


representing the transverse shear strain and 


(3.4.17) 


d A r a 

~dT 


(3.4.18) 
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representing the curvature strain are computed independently from the 
membrane strain A 71 as follows. 

The deformational rotation increments A r a correspond to the rota- 
tional vector parameterization of AS. The computation of (3.4.17) would 
thus require the nonlinear exponential inverse (3.2.13), and the computation 
of (3.4.18) would require an even more complicated derivative of (3.2.13). 
A simpler and more efficient computation can be achieved by incorporating 
Rodrigues parameters to represent A r a and its derivative. The Rodrigues 
parameter representation of the angular variation (3.2.18) is interpreted in- 
crementally as 

Ar a = - 2 ^ 2 [ I + 0 T ] A0 t . (3.4.19) 

The parameters 0 t and A0 ( , which represent the S n+1 and AS matrices, 
respectively, can be easily computed from (3.2.15). In addition, as the S 
matrices are nodal orientations relative to a common elemental reference 
frame T, the corresponding nodal Rodrigues parameters are defined with 
respect to a common basis and thus can be properly interpolated to achieve 
a discrete evaluation of the elemental strains (3.4.17) and (3.4.18). With 
this introduction, the incremental transverse shear and curvature strains 
are computed as follows. 

The elemental shear strain can be computed in two ways. A set 
of Rodrigues parameters can be computed at each of the element nodes as 
the matrices S n+1 and AS are nodal quantities themselves. The algebraic 
average of these nodal parameters, consistent with the standard underinte- 
gration of transverse shear to prevent locking, can then be used to evaluate 
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(3.4.19). A more accurate method interprets the underintegration in a geo- 
metric sense by defining one set of Rodrigues parameters from the relative 
orientation of the mid-element cross-section. This mid-element shear matrix 
S a is defined as 

S a = R a T t 

where the matrix R 0 corresponds to the geometric average of the total 
cross-section orientation matrices at the element nodes and T is the ele- 
ment convected frame. The mid-element total orientation matrix R a can be 
immediately found using the normalized average of the nodal Euler param- 
eters as defined in (3.3.17). The Rodrigues parameters 0 t and A0 t used to 
evaluate (3.4.19) Eire then obtained from the matrices 


S” +1 = R” +1 T n+l3 


As a = s" s: +i 


(3.4.20) 

(3.4.21) 


respectively. 

The elemental curvature (3.4.18) is computed by differentiating 
(3.4.19) as 


dAr a 

dZ 


rref t / + ® r l 


dAQ t 

dZ 


+ 


2 dQ t 

1 + 0 ? dZ 


( 1 + 0? )2 I 1 + ® T ] A0 ‘ 


A 0 t 


(3.4.22) 


The derivative terms Eire approximated over Ein element of length t as 

de t 1 


dZ 

d AQ t 


t 

1 


( ©.. - ©I. ) 


- I ( A®ii - A0,, ) 
using the nodEil Rodrigues parameters 0 t . and A0 t ., i = 1,2. The non- 
derivative terms are evaluated using the mid-element shear matrices defined 
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in (3.4.20) and (3.4.21). This choice is obtained using the same concepts 
introduced in (3.3.17) and (3.3.18) to evaluate a total curvature from the 
Euler parameters of the nodal orientations. 

Thus completes the computational details of the incremental strain 
evaluations. The strain computations are given by (3.4.15) or (3.4.16) for 
the incremental membrane strain, (3.4.17) and (3.4.19) for the incremental 
transverse shear strain, and (3.4.18) and (3.4.22) for the incremental curva- 
ture. These accurate incremental strain computations are a vital ingredient 
of the stress update procedure detailed in Section 2.4. 

3.5 Static Solutions 

The ability of the computational procedure to effectively model large 
deformations is demonstrated as follows. Static equilibrium configurations 
corresponding to given external forces are compared to results in the existing 
literature. A solution to the nonlinear static equations 

S ( u n+1 , R n+1 ) = f? x ? (3-5.1) 

for the configuration which maintains equilibrium between the internal 
stresses and externally applied forces can be obtained from a Newton- 
Rhapson iteration procedure. The Newton-Rhapson technique requires the 
consistent linearization of the internal force as derived in Section 2.5. Given 
the geometric and material tangent stiffness matrices A G and A M repre- 
senting the linearized internal force, the equations 

[* G + * M w, {£},„„ = < 3 - 5 - 2 > 

r n+i = s ( u n+i , R n+i ) - 


(3.5.3) 
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are solved at the k + I th iteration for the incremental translations A u and 
rotations Aa. From this incremental solution, the nodal translations are 
updated via 

“” + <M-i> = “" + W + Au . (3.5.4) 

and the nodal orientations are updated via 

R + (jt+i) = ex P ( ) R • (3.5.5) 

The procedure is initiated with a prediction for ( u n+ Q , R n+ o ). The 
residual is then reevaluated with the updated solution, and the procedure 
continues until the norm of the residual converges within a given range of 
zero. A discussion of the tangent stiffness matrices within (3.5.2) and static 
results obtained from the Newton procedure are given as follows. 

3.5.1 Discrete Tangent Stiffness Matrices 

The finite element method is used to obtain discrete representa- 
tions of the material and geometric tangent stiffness operators presented 
in (2.5.29) and (2.5.31) in the same manner as discussed in Section 3.3 for 
the finite element discretization of the internal force operator. The discrete 
elemental matrices are given as 

8K m = {£ Ul 6u 2 6a ! 6a 2 } T 

for the material tangent stiffness where the partitions K*f. are given as 

‘ T T Ci T -T t C x T' 

. -T T Cx T T t CxT 


TfM _ - 

K " - 1 


jsM I'M 
A 11 12 

k mT i< m 
1X 12 “ 22 



' A ux ' 


Au 2 

i 

Aoi 
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-T t C x ?x 5f -T T C x i x 5f 

k m = - 

2 [ T t Cl »1 5f T t Ci ?! 5f 

iSi ?f Cx ilS r f5i ? f Cx ixS? 

-\SxC 2 5f + i 5! C 2 5f 

sym | 52 if Ci ii5f 

+ 7 52 C 2 Sj 

Ci = Diag ( £M,CA,GA ) 

C2 — Diag ( GJ^El2yEI^ ) 

It can be seen that this material tangent stiffness matrix is simply a trans- 
formation of the standard linear Timoshenko beam stiffness matrix as 

K m = T/ Ktm T k 

where 

T k = Diag ( T,T,5i,5 2 ) 

7 Ci -i Ci Ci ?i Cx 'n 

-i Ci i Ci \ Cx 11 \ Ci 11 

-\i T xCx -\iJCx {if Ci h i'lj Cx 11 

\ if Cl i If Cl i if Cl *1 i if Cl *1 





where 
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0 0 0 

0 0 0 

0 \c 2 -\c 2 

0 \c 2 -\C 2 

As seen from the above decomposition, Ktm possesses the correct number of 
singularities to correspond to the six rigid body modes. Thus, the material 
stiffness matrix K M is not rank deficient when reduced integration methods 
are introduced for the rotational degrees of freedom associated with the 
transverse shear strain. 

Likewise, the geometric element stiffness kernel is given as 
SI\ G = { Sill & u 2 $ a l 2 

where the partitions K G . become 






T t ( \M* - \N* ) Sj T t ( -\M* - \N * ) Sj ' 
T T ( -i M* + |iV* ) Sj" T t ( \ M* + \N* ) S% 


Si J Si -Si J Si ' 

+ 

-S 2 J Si S 2 J Si 
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1 p! ( -0* - o* T ) Sj Si ( -o* + o * T ) s 2 t 

2 |_ S 2 ( O* - 0* T ) Sf S 2 ( o* 4- 0* T ) S 2 r . 

The geometric stiffness is simply 

I\ G = Tj I<tg Tk 

where 

I 7 i M* —j M* 

i J ^ J M* \ M* 

\ m* t m* t f J j 

-i m* t i m* t -4- J 4- J 
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Again, it is seen from the above decomposition that the geometric stiffness 
also contains the correct number of singularities to correspond to the six 
rigid body modes. The geometric stiffness matrix is not rank deficient when 
the reduced integration techniques are employed. 

Two remarks are made in regard to the above derivations. First, it 
is noticed that the elemental curvature k has been neglected. The transfor- 
mation from convected coordinates to body frame coordinates of the virtual 
rotation spatial derivative is effected via (2.3.26) as opposed to (2.3.27). In 
this manner, the spatial derivative employed is as witnessed by the con- 
vected reference frame rather than the cross-section reference frame. For 
static analysis, better convergence results from this interpretation. For dy- 
namic analysis however, the transformation (2.3.27) is retained to coincide 
with the analogous angular velocity representation. As such, the interaction 
between the curvature and angular velocity is accurately represented. Sec- 
ond, the tangent stiffness matrices are symmetric. Previous works achieve 
symmetric tangent operators by incorporating specific rotational parame- 
terizations within the formulation [34-35]. The additional complexity re- 
garding the potential symmetry in these works arises due to the use of the 
unsymmetric first Piola-Kirchoff stress representation, whereas the present 
formulation directly achieves symmetric tangent matrices due to the use of 
the Cauchy stress representation. 

3.5.2 Example Problems 

In this section, two numerical examples are given to demonstrate 
the effectiveness of the formulation in modeling large deformations. The 
results are compared to those reported in the literature. 
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In the first example, a cantilever beam is subjected to a concentrated 
end moment M about the b 3 axis as shown in Figure 3.3a. The beam is 
modeled with twelve linear elements. The total load of M = 2 m is applied 
in ten equal load increments, and an average of 8 iterations is required at 
each increment to reach a converged configuration. The deformed shapes of 
the be am at each load increment are shown in Figure 3.3b. The tip position 
coordinates at each load increment are compared to results reported in [27] 
in Table 3.1 and also the true arc segments of a circle of the appropriate 
radius of curvature. 

In the second example, a cantilever 45-degree bend located in a 
horizontal plane is subjected to a concentrated vertical static end load as 
shown in Figure 3.4a. The bend, of radius 100, is modeled with 8 linear beam 
elements. A total load of 600 is applied in twelve equal load increments, and 
an average of 10 iterations is required at each increment to reach a converged 
configuration. The deformed shapes of the beam at each load increment are 
shown in Figure 3.4b. The tip position coordinates at each load increment 
are compared to results reported in [24,32] in Table 3.2. 

The results are given as follows. 
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Figure 3.4a Cantilever beam subjected to end moment: problem data. 
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Figure 3.5a Cantilever 45-degree bend subjected to end load: problem data. 
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3.6 Concluding Remarks 

The finite element representation of the equations of motion derived 
in Chapter II has been presented. Much attention was focused on the dis- 
crete form of the internal force, an accurate computation of which remains 
key to any static and/or dynamic simulations. The numerical procedure for 
the computation of the nonlinear internal force, as developed from the con- 
tinuum formulation, was proven to be invariant to arbitrary rigid motions of 
the beam while effectively modeling the strains due to stretching, transverse 
shear, torsion, and bending deformations. Particular developments include 
the discrete computation of rotational matrices embedded within the formu- 
lation, a discrete computation of the element curvature based on the Euler 
parameter functional relation, and incremental strain computations to effect 
the stress update procedure. The fundamental property motivating these 
developments is the condition that quantities be computed from the nodal 
degrees of freedom in a manner that is invariant to any arbitrary rigid body 
motions implicitly contained within these variables. The success of the in- 
cremental strain computational procedure hinges on a proper extraction of 
the rotation increments from the various rotational matrices embedded in 
the formulation. 

The developments of multibody dynamic solution procedures is pre- 
sented in the following chapter. The present internal force computation is 
interfaced with these developments to achieve effective dynamic simulations. 



CHAPTER IV 


MULTIBODY DYNAMICS SOLUTION PROCEDURES 
4.1 Introduction 

The formulation modeling the large rotation dynamics of a flexi- 
ble be am has been developed in Chapter II. The algorithmic treatment of 
the beam internal force such that the strain computations are preserved 
under the large rotations was detailed in Chapter III. The purpose of the 
present chapter to extend these developments to the framework of flexible 
multibody systems. To achieve an effective simulation methodology, an au- 
tomatic formulation of the equations of motion for multibody systems is 
necessary. In addition, efficient numerical solution procedures need to be 
developed for the solution of these equations. The numerical procedures for 
multibody systems must include a proper treatment of the spatial rotations 
and the constraint conditions. The analytical and computational formalisms 
presented in Chapters II and III are highly conducive to these goals. The 
present chapter develops an effective numerical solution strategy for flexible 
multibody dynamic systems. 

To formulate the equations of motion for multibody systems, two 
basic approaches have been discussed in the literature. One approach an- 
alytically eliminates the constraint conditions a priori to achieve a mini- 
mal set of second order differential equations [69-76]. In these recursive 
type formulations, relative coordinates between bodies are used to define 
the kinematic relations, and modal coordinates are used to define the local 
deformations. However, when the kinematic constraint equations are intro- 
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duced in the early stages of analysis, the development of the equations of 
motion becomes unnecessarily complicated. Equations with a high degree 
of nonlinearity and low degree of sparsity result when the relative coordi- 
nates are employed. In contrast, a second approach introduces absolute or 
Cartesian coordinates to model the kinematic relations of each individual 
body in the system. The Lagrange multiplier technique is then employed 
to couple the algebraic constraint equations with the differential equations 
of motion; the result is a maximum set of differential and algebraic equa- 
tions. This approach has been applied to the floating frame formulations of 
beam dynamics where the displacements of elastic bodies are modeled by 
superposing local elastic deformation coordinates onto rigid-body reference 
coordinates [8-11,16,77]. To employ the approach within multibody appli- 
cations, the assumed elastic displacement field must be consistent with the 
nonlinear constraint equations. As simple free-body modes are not suit- 
able to describe the deformation shape of many constrained components, 
consistent displacement fields are defined via a proper choice of the body 
reference frame [8]. To better account for reactions due to joint connections, 
static correction modes are introduced to model the elastic displacements of 
flexible components [10-11]. The important issue thus becomes the proper 
choice of modal coordinates to represent constrained components. 

The present beam formulation alleviates this issue and is extended 
into the field of multibody dynamics as follows. As discussed in Chap- 
ter II, the beam components are modeled using inertially-based degrees of 
freedom which embody both the rigid and deformation motions. As these 
degrees of freedom are kinematically of the same sense as the physical coor- 
dinates of rigid body components, the equations of motion for tin arbitrary 
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configuration of flexible beams and rigid bodies can automatically be gen- 
erated in terms of an identical set of physical coordinates. Incorporation 
of multibody system constraints via the Lagrange multiplier technique be- 
comes straightforward. A salient feature of this type of formulation is that 
numerical solution procedures for treatment of multibody systems can be 
directly applied to the generalized coordinates of both the rigid and flexi- 
ble components. In comparison to the floating frame types of formulations, 
concerns of correct modal representation of the deformation are eliminated 
as well as the separate numerical treatment of the contrasting characters of 
the reference and elastic coordinates. 

Given the flexible multibody system model, specialized numerical 
procedures must be developed for the solution of the equations of motion. 
As the equations of these articulated systems inherently involve large rota- 
tions and constraint conditions, numerical procedures must include treat- 
ment of these effects as well as the integration of the generalized coordinates. 
In the context of present type flexible beam dynamics formulations, the ex- 
tension of the implicit Newmark algorithm to include proper treatment of 
the nonlinear rotation fields was proposed in [32], and consequently used in 
the works of [34,36]. To treat the spatial rotations, the equations of motion 
axe written in terms of a paxticulax rotational parameterization and then 
discretized by an implicit algorithm. The configuration orientation and an- 
gular velocity are updated from the solution for the rotational parameters by 
implementing discrete counterparts to the functional relations overviewed in 
Section 3.2. In contrast to these implicit integration techniques, the present 
work develops explicit integration techniques for use within the rotational 
dynamics context [44]. For this purpose, a two-stage modification of the 
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central difference algorithm is presented to integrate the translational co- 
ordinates and the angular velocity vector. Once the angular velocities axe 
obtained, the angular orientations are updated by incorporating an implicit 
trapezoidal algorithm to solve the kinematical relation in terms of the Euler 
parameters. 

In addition to the integration of the generalized coordinates and 
the update of the angular orientation, solution techniques must be devel- 
oped to enforce the algebraic constraint conditions for multibody equations 
of motions formulated by incorporating Lagrange multipliers. For this pur- 
pose, techniques have been developed which eliminate excess or dependent 
variables numerically by incorporating the singular value decomposition [78- 
80], the coordinate partioning technique [81], the zero eigenvalue theorem 
[82], and other matrix procedures [83-85]. Another approach applies an im- 
plicit integration algorithm to the generalized coordinates of the constraint- 
augmented equations, and then solves for both the generalized coordinates 
and the Lagrange multipliers simultaneously [86-88]. However, in some in- 
stances, numerical difficulties were encountered from this discretization as 
ill-conditioned matrices and artificially stiff problems resulted. Another si- 
multaneous solution procedure attempted to alleviate this problem by con- 
structing a modified differential equation to model the constraint equations 
[89-90]. Drawbacks to the procedure included a somewhat ad-hoc selec- 
tion of certain factors that influence the accuracy of solutions and a lack 
of a positive error control on the constraint violations. To address these 
drawbacks, a technique which implicitly integrates an alternative stabilized 
companion differential equation for the constraint forces was proposed in 
[91-92]. A principal feature of this method is that the errors committed in 
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each constraint condition decay with a characteristic time scale associated 
with the constraint force. In addition, the method is conducive to an at- 
tractive modular software implementation as the solution of the constraint 
forces can be performed separate from the integrations of the generalized 
coordinates. This stabilized constraint solution procedure is employed in 
the present work. 

The multibody dynamics solution procedure adopted for the flexible 
multibody system equations of motion is presented in detail in this chap- 
ter. Section 4.2 will briefly describe the extension of the beam formulation 
into the equations of motion for multibody systems. Section 4.3 will dis- 
cuss the application of integration procedures, both implicit and explicit, 
in the context of finite rotation dynamics. To complete the multibody dy- 
namic solution, Section 4.4 will present a unique procedure to update the 
configuration orientation from the angular velocity solution and Section 4.5 
will discuss a separate computational procedure for the Lagrange multipli- 
ers. Finally, example problems demonstrating the software capabilities are 
presented in Section 4.6, and conclusions are given in Section 4.7. 


4.2 Multibodv System Equations 

The present formulation of spatial beam dynamics can readily be 
incorporated into a general multibody dynamics methodology. The discrete 
equations of motion of an unconstrained flexible beam as derived in (3.3.12) 


are rewritten below as 
m 0 ( it \ 

0 J I u I 


wJw 


(3.3.12) 


where the subscript d representing nodal discretizations has been dropped. 
It is recalled the beam nodal degrees of freedom (u,u) corresponding to 
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the inertial coordinates of the beam neutral axis position and the body- 
fixed coordinates of the cross-sectional angular velocity. These degrees of 
freedom are of the same sense as those of a rigid body, namely the inertially- 
based translational position of the center of mass and the angular velocity 
of a moving reference frame orienting the body. Thus the equations of 
motion (3.3.12) can equally well represent a rigid body system by setting the 
internal force S equal to zero. The unconstrained equations of an arbitrary 
configuration of flexible beams and rigid bodies can be written in terms of 
one set of kinematical coordinates denoting both the nodal coordinates of 
the flexible members and the physical coordinates of the rigid bodies. 

It remains to augment the nonlinear algebraic equations which 
model the contact conditions of the various bodies of the multibody system. 
The Lagrange multiplier technique augments these constraint equations to 
the unconstrained differential equations of motion. Two types of constraint 
conditions exist, holonomic or configuration constraints and nonholonomic 
or motion constraints. A holonomic constraint condition can be expressed as 
a definite relation between the displacement coordinates, whereas a nonholo- 
nomic constraint condition can only be expressed as a relation between the 
differentials of the coordinates and not as a finite relation between the coor- 
dinates themselves [98]. A set of algebraic equations representing holonomic 
constraint conditions between the displacement coordinates u are written as 

= 0 (4.2.1) 

and the differential of these constraints is given as 

d<b H 

S$ H = — — 6u = B h Su = 0 . (4.2.2) 

OU v ' 
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A set of non-integrable equations concerning nonholonomic constraint con- 
ditions between the virtual displacements and rotations 8u,8a axe written 
as 


8$ n ( u,8u,R,8a ) 



0 


(4.2.3) 


Given these constraint conditions, the virtual work expression (2.3.2) for the 
unconstrained system is modified by the inclusion of the virtual work which 
enforces the constraints through the Lagrange multiplier technique as [98] 


8F i + 8F S + A h • 8Q h + \n ■ 6*n = 8F^ + 8F ^ 


By incorporating the Lagrange multipliers A, the virtual displacements and 
rotations of the generalized coordinates can be treated as independent vari- 
ations. The equations of motion for constrained flexible multibody systems 
are thus obtained as 


m 0 
0 J 


]{*} 


+ B r A 


■ (£} ■ 


(4.2.4) 


In the above equation, the notation contains both the holonomic and non- 
holonomic constraints in the constraint force vector A as 


B = 


Bh 

Bn 


_ / 

“ \A N 


(4.2.5) 


and the right-hand side vector contains the remaining force-type terms as 


{£} - 


ft _ 5 e 

f b — S b — CjJuj 


} 


(4.2.6) 


To achieve a determined system of equations, the equations (4.2.4) must be 
solved in conjunction with the constraint equations. 

The B matrix in (4.2.4), termed the constraint Jacobian matrix, is 
deduced from the kinematic relationship between the bodies of the system. 
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To model these kinematic relations, conditions can be directly imposed on 
the nodal degrees of freedom of two separate flexible beams, or between 
a beam nodal degree of freedom and a rigid body modeling a portion of 
a more complex joint. Typical Jacobian matrices modeling standard joints 
used in multibody systems, such as universal, spherical, revolute, and trans- 
lational joints, have been derived in [99]. Given the Jacobian matrices for 
the joints and the beam connections, the models representing an arbitrary 
assemblage of articulated flexible and rigid components can be constructed 
in a systematic manner. 

4.3 Integration of Generalized Coordinates 

The solution of the generalized coordinates via the application of 
numerical integration techniques to the multibody equations of motion are 
discussed next. To present the application of integration methods to prob- 
lems concerning nonlinear rotational dynamics, the equations of motion at 
time t n+1 are written in terms of a generalized coordinate d which represents 
both the displacement coordinates u and a suitable rotational parameteri- 
zation as 

M (d n+ \d n+1 ,d n+1 ) + S ( d n+1 ) = /( t n+1 ) (4.3.1) 

where d n+1 = d ( i n+1 ). In the above equation, M represents both the 
translational and rotational terms of the inertia operator as 

. (4 - 3 - 2) 

S represents the nonlinear interned force, and / represents the external force 
including any constraint forces / = f ext — B r A. For the present analysis, 
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the constraint forces are assumed to be given; a separate solution technique 
in which the constraint forces may be obtained in conjunction with the 
solution of the generalized coordinates is discussed in Section 4.5. Numerical 
algorithms axe implemented to advance the generalized coordinate solution 
to a discrete point in time t”'*' 1 = t n + h given that the solution at the 
past step in time t n has been found. 

To solve equations of motion involving nonlinear rotational dynam- 
ics, implicit integration techniques were first considered [32,34,36]. The 
advantage of the implicit algorithms is the exceptional stability properties, 
while the disadvantages include the implementation complexities inherently 
associated with nonlinear problems. Most of the implicit methods that are 
commonly used, i.e. the trapezoidal rule, Newmaxk, Houbolt, Wilson, and 
others, are unconditionally stable for linear problems. However, to achieve 
this unconditional stability, the solution of a simultaneous system of equa- 
tions is required at each discrete time step where the dynamic solution is 
sought. To apply implicit techniques to nonlinear problems, a linearization 
of the equations with respect to the generalized displacement coordinates 
is necessary. For the linearization within the context of finite rotations, 
the rotational degrees of freedom are parameterized by the rotational vec- 
tor defined in Section 3.2 [32]. Implicit integration methods are then used 
to obtain the rotational vector and its time derivative at the discrete time 
stations, after which the angular orientation and the angular velocity are 
updated accordingly. 

The alternative to the implicit formulas is the explicit formulas. 
Explicit formulas are attractive from an implementation standpoint as the 
solution can be achieved from simple vector operations. However, explicit 
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algorithms retain stability only within a restricted step size between dis- 
crete time stations. To use explicit methods within the context of finite 
rotations, the angular velocity vector is first deduced via an integration of 
the angular acceleration [44]. A suitable discretization of the Euler param- 
eter representation of angular velocity is then used to update the rotational 
orientation. 

To demonstrate the above concepts, the application of the implicit 
midpoint rule is discussed next in Section 4.3.1, while the explicit techniques 
axe discussed in Section 4.3.2. 

4.3.1 Implicit Integration Techniques 

A traditional algorithm to solve (4.3.1), termed the Newmark algo- 
rithm, is given as 

d n+1 = d n + hd n + ( i - 0 ) h 2 d n + 0 h 2 d n+1 

d n+1 = d n + ( 1 - 7 ) h d n+1 

where 0 and 7 are free parameters whose choice results from a compromise 
between computational simplicity, stability, and accuracy requirements. The 
choice 7 = | , 0 = j yields the most accurate integration scheme 

with unconditional stability. The resulting algorithm, termed the average 
constant acceleration scheme or trapezoidal rule, can be rewritten as 

d n+i = d n + £ ( d n + d n+1 ) 

£ 

d n+1 = d n + £ ( d n + d n+1 ) . 

A variation of this scheme, termed the midpoint rule, is given as 


d n+l 

= d n + hd n ^ 

(4.3.3a) 

d n+1 

= d n + hd n+ * 

(4.3.36) 
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where a computation of the velocities and accelerations at a half time station 
are evaluated from the application of 

d n+ ? = d n + ^ d n+ * (4.3.4a) 

d n +i = d n + - d n+ * (4.3.46) 

2 

to the equations of motion. By evaluating the velocities and accelerations 
at the half time station for use in the full station displacement and velocity 
computation as opposed to using a numerical average as seen in the trape- 
zoidal rule, the midpoint algorithm yields a more accurate computational 


scheme. 

The midpoint solution is obtained by using the Newton-Rhapson 
iteration procedure to solve the nonlinear residual equations 
r n+ * = M ( J n+ M n+ M n+ * ) + S(d n+ M - / n+i 

(4.3.5) 

= 0 


To this end, a linearized set of equations 

n+ ^ 


E t 


(fc+D Ad(*+i) 




(fc+D 


E 


d r n+ ^ 
dd n+ ? 


( 4 . 3 . 6 ) 


are solved at the k + I th iteration for the incremental displacements Ad 
which are then used to update the solutions obtained at the k th iteration 
via a discrete evaluation of 


( 4 - 3 - 7 > 

until a convergence criterion is reached. The iterative procedure is started 
with a prediction for usually taken to be d n . The configuration 

update (4.3.7) symbolically designates the vector update 


u 


n+l 

(*+l) 


U 


n+l 

(*) 


+ A u 


(4.3.S) 
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for the translational displacements u and the rotational orientation update 

R " + (‘* + .> = «J> ( ) R ” + (‘ 4) (4-3.9) 

for the rotational vector parameterization 0. The system angular velocity 
must be updated from the rotational vector and its derivative in an appro- 
priate manner via 


w n+1 = T ( 0 (n+1) ) 0 n+1 


where the operator T ( 0 ) has been defined in 3.2.17. 

To obtain the linear equations, the incremental solution matrix E 
is derived from (4.3.5) and (4.3.6) as 

9M n+ i d d n+ * 3M n+ i d d n+ * 


E = 


d d n+ 5 d d n+ 7 


+ dd n +? dd n+ s + 
d M n +* d 5 n+ 2 


where the terms 

d d n+ i 
d d n+ ? 


1 _ 

8 2 


d d^2 
d d n+ i 


d d n + j 


1 

8 


+ 


d d n+ % 


6 - 5 


(4.3.10) 


(4.3.11) 


axe identified from the integration algorithm (4.3.3). The partied derivatives 
of M and S represent coefficients of the first-order Taylor-series expansion 
of the respective operators about a past equilibrium configuration at time 
t n as 


M ~ 
5 ~ 


M n + 

S n + 


d M 

d d 
d S 


Ad + 


d d 


Ad 


d M 

~dl 


Ad 


+ 


d M 
~d~d 


Ad (4.3.12) 


(4.3.13) 


where Ad = d — d n . Derivation of the linear expansions in the context 
where d = ( tt,0 ) requires techniques consistent with the particular 
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translational and rotational update procedures of (4.3.8) and (4.3.9); the 
derivations for the present beam formulation were presented in Section 2.5. 
The consistent linearization of the internal force as derived in Section 2.5 
and Section 3.5.1 is given as 

5 ~ S n + [ K g + K M ] { ^0 } (4.3.14) 

where the geometric and material tangent stiffness matrices I\ G and K M are 
deduced from a finite element discretization of (2.5.28) and (2.5.30). Like- 
wise, the consistent linearization of the translational and rotational inertia 
as derived in Section 2.5 is given as 


M m ~ M” + m Att , (4.3.15) 

Mj ~ M] + J A0 + C G A0 + K c A0 (4.3.16) 

where C G and K c represent discrete versions of the gyroscopic damping 
and centrifugal stiffening matrices given by (2.5.14) and (2.5.15). From the 
above linearizations, the incremental solution matrix is seen to be 

1 \m 0] 1 f 0 0 1 [0 0 1 k g , k m 

E = P [o j\ + « [o H + [o I< c \ + A + A 

This solution matrix E is configuration dependent and thus must 
be reformulated at least at the beginning of each time step for a modified 
Newton technique, and at each iteration if a full Newton technique is to be 
used. In addition, the updates of the configuration orientation and angular 
velocity from the rotational vector parameterization require additional cum- 
bersome computations. It is seen that the price to pay for the unconditional 
stability of the implicit algorithm is the complexity of the computational 


scheme. For these reasons, explicit integration procedures have been chosen 
for the present study. 


4.3.2 


where 


Explicit Integration Techniques 

The central difference explicit integration algorithm is written as 



= d"”* + hd n 

(4.3.17a) 

d n+i 

= d n + hd n+ * 

(4.3.176) 

d n+1 

= M- 1 Q ( d n+1 , d n+1 ) 

(4.3.17c) 

{?} 

- (?) - Ul 

(4.3.18) 


In this algorithm, the displacements ( d ) are obtained at the full time steps, 
( • • ■ , n, n -f- 1, n + 2, • • • ), while the velocities ( d ) are obtained at the half 
time steps, ( •••,« — j,n + |,n+|,-*- ). Given the integrated displacements 
and velocities, the accelerations are obtained at the full time steps from the 
equations of motion. The solutions are thus advanced to the next time step 
using an extrapolation involving only the previous solution vectors, and the 
implementation of explicit integration techniques to the nonlinear equations 
of motion is straightforward. 

To effect the explicit formulas within the context of finite rotations, 
the algorithm (4.3.17a) is used to obtain a solution for the angular velocity. 
The rotational orientation can then be deduced from the angular velocity 
via a separate numerical procedure discussed in Section 4.4. However, com- 
putational difficulties arise when (4.3.17a) with ( d => u> ) is employed to 
solve the rotational equations 
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typical of multibody dynamics. From these equations, it is seen that in order 
to compute w n , it is necessary to have w n . However, due to the inherent 
nature of the algorithm, only u j n ~2 is available. To alleviate this problem, 
the following approximations may be introduced: 

u n ~ u n ~* (4.3.20) 

w n ~ - ( u n+ > + ) . (4.3.21) 

2 

It will shown that the first approximation results in a computationally un- 
stable algorithm, while the second approximation remains stable but does 
not retain an explicit nature of computations. These difficulties motivated 
the development of a modified central difference algorithm termed the two- 
stage staggered procedure [44]. This algorithm is based on an interlaced 
application of the central difference algorithm such that the displacements 
and velocities are advanced one-half time step at a time. As standard of 
explicit integration methods, the new algorithm remains stable under a re- 
stricted step size and is based on vector operations. The stability analysis 
of the algorithms follows. 

4.3.2. 1 Standard Central Difference Implementation 

For a linearized stability analysis of the central difference algorithm 
(4.3.17) with approximations (4.3.20) and (4.3.21) as applied to rotational 
equations typical of multibody dynamics, the linearized equations of dy- 
namics 

J A0 + C G A0 + ( K m + K g + K c ) A0 = f u - Mj + Mj 

as derived in Section 2.5.1 will be examined, the torque- free motion of a 
single spherical rigid body in which the principal inertia terms are equal, 
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the above linearized equations reduce to 


I 0 + u 0 = 0 . (4.3.22) 

This simple case is used to demonstrate important stability characteristics 
as follows. The effects of the geometric, material, and centrifugal stiffnesses 
K G , K M , and K c are the deduced from numerical experiments to asses the 
numerical stability of flexible systems. 

To first examine the effects of the combination (4.3.17) with (4.3.20) 
on the linearized equations (4.3.22), we apply (4.3.17) with ( d =$> 0 ) to 
(4.3.22) to yield 


0"+i = _ huQ n (4.3.23) 

0" +1 = 0” + h 0” + * . (4.3.24) 

The approximation (4.3.20) with ( u> =» 0 ) is then introduced to give 

©"+* = [/-fcwje—i (4.3.25) 

0 n+1 = 0" + 6 0 n+ 5 . (4.3.26) 


To assess the stability characteristics, these equations are rewritten in the 
form of a difference equation in terms of 0 as 


( 0 n+1 - 20 n + 0"" 1 ) + h u ( 0 n - 0 n_1 ) = 0 . (4.3.27) 

Computational stability is determined by seeking a nontrivial solution of 
the form 

0 n+1 = sQ n = s 2 0 n_1 . (4.3.28) 

For a stable algorithm, it is necessary that 


I a | < 1 


(4.3.29) 
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From (4.3.28), the characteristic equation of (4.3.27) becomes 

det | ( 6 - 1 ) 2 + huj(s - 1)| = 0 . (4.3.30) 

In order to test the stability requirement (4.3.29) on the above equation, 
the unit circle (4.3.29) is transformed into the entire left-half plane of the 
z-domain via 

5 = J-t-f (4.3.31) 

1 — z 

such that (4.3.29) is equivalent to 

Re(z ) < 0 . (4.3.32) 

The characteristic equation (4.3.30) is evaluated by introducing (4.3.31), 
and the z-polynomial equation 

h 2 ( u> 2 + a >2 + wf ) z 2 — 2 h 2 ( u 2 + + u/f ) z 

(4.3.33) 

4- [ 4 + h 2 ( ) ] = 0 

results. This equation must be examined for possible unstable roots. For 
this purpose, the Routh-Hurwitz stability criterion. [95] is employed. In 
order for all of the roots of a second-order polynomial to lie in the stable 
left-hand plane, each of the coefficients of the polynomial must be non- 
negative. It is immediately seen that the coefficient of the z-term in the 
above is always negative. Thus, for the present rigid body example, the 
central difference algorithm becomes an unstable integration method when 
the conventional velocity approximation is incorporated. This property is 
witnessed in numerical simulations of flexible multibody systems when the 
given procedure is used to integrate the equations of motion. For this reason, 
this procedure is not suitable for multibody dynamics applications of the 
form (4.3.19). 
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The alternative velocity approximation (4.3.21) as implemented into 
the central difference algorithm (4.3.17) corrects this instability problem. In 
analyzing this procedure, the substitution of (4.3.21) with ( u =>• 0 ) into 
(4.3.23) gives the counterparts to (4.3.25-4.3.26) as 

( i + | a ) e"+i = (/-£«) e-* 

0 n+1 = 0 n + h 0 n+ ? 

The equivalent difference equation for the above algorithm is 

( 0" +1 - 20 n + 0”’ 1 ) + £ u ( 0 n+1 - 0 n_1 ) = 0 . (4.3.34) 
The characteristic equation of the above 

det | ( s - 1 ) 2 + ^ u ( s 2 - 1 ) | = 0 (4.3.35) 

is again transformed to the z-domain via (4.3.31). In this case, the z- 
polynomial equation becomes 

.z 4 [ 16 z 2 + h? ( + cc »2 + ^3 ) ] = 0 (4.3.36) 

It is easily seen that for this case roots are either zero or purely imaginary. As 
such, the central difference algorithm combined with the alternative velocity 
approximation is unconditionally stable for the simple rigid body example. 
For systems of flexible bodies, the stability limit of the above algorithm is 
governed by the maximum elastic frequency of the beam component which 
can be determined from the eigenvalues of the tangent stiffness matrices 
K G and K M or by the centrifugal stiffness of the tangent mass K c . From 
computational experiments, it has been seen that the stability limit closely 
corresponds to that as predicted from the elastic frequencies of linear Tim- 
oshenko beam stiffness matrices. 
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The one drawback to the use of this last algorithm is that implicit 
computations are required to determine the velocity when the equations of 
motion contain velocity dependent force terms other than diagonal damping. 
To show this, the algorithm (4.3.17a) with (d =* w ) and the approximation 
(4.3.21) are applied to the nonlinear equations of motion (4.3.19). The 
following nonlinear residual equation 


r 


n+i 


J + w""* ) J ( u; n+ * + u; n ~* ) 

4 

- J - hf u = 0 

(4.3.37) 


results. The Newton-Rhapson technique must be employed to solve the 
above equation for w n+ 2 . For this case, the following iterative sequence 


E k Aw = ~r n+ l (4.3.38a) 

w n+ L +1 = u n k + Aw (4.3.386) 


results, and the solution matrix E is easily derived as 


E 


J\ | {J 3 ~ J 2 ) ^3 4 (^3 — J 2 ) &2 

T (Jl — J 3 ) 0)3 J 2 4 {J\ — ^3) 

4 (J 2 — <^l) £ (J 2 ~ <A) Wi J 3 

w = ( w""*”* + u) n * ) 


The implicit computation is only necessary to obtain the angular velocities of 
multibody system equations. The translational displacements and velocities 
are found from the standard application of the central difference algorithm, 
and the rotational orientation is found from the angular velocities with a 
procedure to be discussed in Section 4.4. The necessity for the Newton- 
Rhapson solution of the angular velocity can be eliminated by the following 
algorithm. 
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4.3. 2. 2 Modified Central Difference Implementation 

An alternative algorithm has been derived to alleviate the problems 
or complications resulting from the velocity approximations which must be 
introduced when the central difference algorithm is used to integrate the 
equations of motion for multibody systems. The troublesome velocity ap- 
proximations can be eliminated by staggering the standard central difference 
update such that the displacements and velocities are advanced one-half time 
step at a time. The algorithm advances the solution to the time station t n+ 2 
given the solutions of the two preceding time stations t n ~ 2 and t n as 


d n +i = d n ~* 

+ 

h d n 

(4.3.39a) 

d n+ * = d n_ 5 

+ 

h d n . 

(4.3.396) 

The algorithm is initiated for time t* given 

initial conditions for time t° as 


follows: 

d* = d° + 

d 1 = d° + hd> 

di = \(d? + d') . 

The result is an explicit computation which remains stable under a restricted 
stepsize. 

The stabilizing effect of this algorithm is demonstrated with the 
same procedure used in assessing the other algorithms. Application of 
(4.3.39) to the linearized equations (4.3.22) yields the difference equation 

( 0 n+1 - 20"" 1 + 0"“ 3 ) + hu(G n - 0"~ 2 ) = 0 

where h = 2 h. The associated characteristic equation is seen to be 

det | ( s 2 — 1 ) 2 + h u s ( s 2 - 1)| = 0 
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from which the 2 -domain stability polynomial is derived as 

h 2 ( + u 2 + ) z 4 + [ 16 - 2 h 2 ( u>j + ) ] z 2 + 

h 2 ( U 2 + U 2 + ) = 0 . 

For a fourth-order polynomial of the form 

az 4 -|-6z 2 -|-c = 0 

to guarantee stable roots, the conditions 

a > 0 , b > 0 , c > 0 , b 2 — 4ac > 0 

must all be satisfied. This is accomplished by restricting the stepsize to be 

\f w 2 + o>2 + ^3 

for the specialized example problem. This linearized stability analysis con- 
firms that the present two-stage explicit procedure in fact remains stable in 
the presence of gyroscopic damping. The stable stepsize for flexible multi- 
body applications is dictated by the above and the standard central differ- 
ence limit 
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where ut is the highest elastic frequency of the beam component. As verified 
by computational experiments, the flexible stability limit can be approxi- 
mated from the linear Timoshenko beam stiffness matrices. 

This algorithm, termed the two-stage explicit procedure, can be ap- 
plied to the nonlinear equations of motion in a straightforward manner; 
(4.3.39a) with ( d => u , u> ) is used to integrate the translational and 
angular velocities and (4.3.39b) with ( d => u ) is used to integrate the 
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translational displacements. The algorithm advances the solution a half 
time station at a time, and requires the evaluation of the right-hand side 
terms Q twice each time step. Included in this evaluation is the computa- 
tion of the internal force as discussed in Chapter III. The trade-off between 
this algorithm and the standard central difference algorithm accompanied 
with the velocity averaging, aside from the step-size restriction, is then a 
consideration between the extra force evaluation per time step required of 
the former versus the simultaneous equation solver required of the latter for 
the iterative angular velocity solution procedure. 

Thus concludes the discussion of the integration of the equations 
of motion for multibody dynamics applications by the central difference 
algorithm. The two stable algorithms which have been discussed, namely the 
standard central difference algorithm including the velocity averaging and 
the two-stage staggered explicit algorithm, were introduced for integration 
of the translational displacements and velocities and the angular velocities 
of multibody systems. As the rotational orientation parameters are not 
directly integrable from the angular velocity vector, a procedure must be 
developed to update the configuration orientation given the angular velocity. 
This is discussed next. 

4.4 Rotational Orientation Update Procedure 

The effective treatment of rotational quantities inherently involved 
in the multibody system equations requires special consideration. As ad- 
dressed in Section 4.3.1, implicit algorithms are based on an integration of 
the displacements and a particular rotational parameterization as the pri- 
mary variables. Given the rotational parameterization, usually taken to 
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be the rotational vector, the configuration orientation is directly updated 
via (4.3.9) or a suitable counterpart. Computationally, this can be accom- 
plished using an exact functional expression for the matrix exponential [32], 
using the Euler parameter representation of rotation matrices after convert- 
ing the incremental rotational vector to Euler parameters [32], or using a 
second-order approximation to the matrix exponential [50]. For the explicit 
techniques, the rotational orientation must be deduced from the angular 
velocity solution. 

A classical rotational update procedure is based on a discrete ap- 
proximation 


R n "K ^ eX p ( h ) R n 


u 


— Cb ( t n + /? h ) 


f n+l 


to the true solution 

R n+1 = exp ( r"fi(r)dr)R- 
of the generating differential equation 


(4.4.1) 

(4.4.2) 


R = u T R . (4.4.3) 

The evaluation of the angular velocity at time t n+5 with /? = | in (4.4.2) 
was determined to provide the highest accuracy [96]. A related result ap- 
plied the generalized midpoint rule to the generating kinematic differential 
equation to yield [97] 

R n +1 ~ ARR n 

AR = {I - /3u* h}- 1 [I + (l - 0)u* h, ] . 
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In order for AR to remain an orthogonal rotation matrix, it was shown that 
for the above expression the value /? = \ is needed. The rotational matrix 
approximation can be inverted as 

0 x-l / . 0 X . 1 


(J- f )-‘(x+ =) = / +r 


[ e + ; ® 2 ] 


+ ( ||0|| 2 /4 ) 1 ’ ' 2 

which is equivalent to the second-order Pade approximation of the matrix 
exponential exp ( 0 ) . Therefore, a rotational update is obtained from a 
true or approximate matrix exponentiation of 0 = u ; n+ 2 h. 

An alternative procedure introduced here is based on the numerical 
integration of the following Euler parameter representation of the generating 
equation 


q = A(u>) q , 
in which A is given by 


A <"> - 2 


-{;} 


(4.4.4) 


0 —IxT 

U} up 


for the body-fixed angular velocity coordinate definitions. The configura- 
tion orientation is obtained from a numerical time discretization of (4.4.4). 
Among several possibilities, the approximation that satisfies the inherent 
constraint condition 

qo qo + q T q = 0 (4.4.5) 

in the discrete sense is the following trapezoidal formula 

( j n + 1 — a n 1 1 

q - r - q - = A(u," + i) - ( + ,") . 

Due to the structure of A, the solution matrix can be analytically inverted 
such that the discrete orientation update is accomplished according to 


?" +1 = i [/ + | A( u " + i)] [/ + | A(w” + i)] 


(4.4.6) 
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where 

D = 1 + ( w\ + w\ + w\ ) . 

The final result is normalized to satisfy the Euler parameter constraint con- 
dition. 

The resulting update (4.4.6) involves only explicit computations and 
is readily incorporated into the two-stage explicit integration procedure. 
This Euler parameter-based solution procedure also interfaces with the in- 
ternal force computational procedure presented in Chapter III. The concept 
of an Euler parameter-based orientation update procedure was introduced in 
[44] for rigid multibody dynamic systems. For the rigid cases, many implicit 
integration algorithms as applied to (4.4.4) resulted in accurate numerical 
simulations. For the extension to flexible systems, it was determined that 
the particular algorithm (4.4.6) must be implemented to achieve a stable 
simulation. This particular algorithm implicitly satisfies the constraint con- 
dition (4.4.5), whereas alternative algorithms do not possess this property. 
Thus, successful simulations of flexible systems are not necessarily immedi- 
ately obtained from existing rigid-body dynamic solution procedures. The 
present formulation, involving the internal force computation of Chapter 
III, the generalized coordinate integrator of Section 4. 3. 2. 2, and the orienta- 
tion update procedure of the present section, leads to effective simulations 
as demonstrated in Section 4.6. The next section discusses the solution 
of the Lagrange multipliers which have been incorporated to augment the 
constraint conditions to the differential equations of motion. 
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4.5 Constraint Force Solution Procedure 

A partitioned solution procedure has been employed to solve the 
the Lagrange multipliers and the corresponding constraint forces separately 
from the solution for the generalized coordinates. A stabilized constraint 
force solution procedure originally developed for rigid multibody systems 
in [91,92] was successfully extended to models incorporating the present 
beam formulation. Again, the success of various existing constraint force 
solution procedures developed from analyzing multirigid body systems is 
not necessarily guaranteed for flexible multibody systems. For this reason, 
a particular implementation of the basic formalism of [91,92] was determined 
to yield effective simulations of flexible multibody systems. 

To effect a partitioned solution of the constraints, a stabilized com- 
panion differential equation for the constraint forces is formed by adopting 
the penalty procedure. The penalty procedure uses the equations 

Xh = ~ $H » X N = - $ N , e — ► 0 (4.5.1) 

e e 

as the basic constraint conditions instead of (4.2.1) and (4.2.3), respectively. 
From (4.2.2) and (4.2.3), both penalty equations can be written in the gen- 
eral form, 



where the notation includes both holonomic and nonholonomic constraints 
with (4.2.5). 

The numerical solution to the companion differential equation (4.5.2) 
is obtained as follows. The constrained equations of motion (4.2.4) are in- 
tegrated once using a general implicit integration algorithm given as 

d n+1 = S d n+1 + g n 
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Two specific algorithms have been shown to work within the context of 
flexible multibody systems; they are the trapezoidal rule ( g n = d n + 
Sd n , 8 = j ) and the forward Euler formula ( g n = d n , S = h ) . The 
result of the integration 

d n+1 = 6 M~ l ( Q n+1 - B t A n+1 ) + g n (4.5.3) 

is substituted into (4.5.2), and a stabilized differential equation 

eA n+1 + SB M~ l B t A n+1 = SB M _1 Q n+1 + B g n (4.5.4) 

is obtained for the Lagrange multipliers. To obtain a discrete update of A, 
the same integration rule is applied to this equation as was used for the 
generalized velocity. The trapezoidal rule 

A n+1 = A" + S ( A" + A n+1 ) 


results in a constraint force update of the form 


(e I + S 2 BM _1 B t ) A n+1 = (el - S 2 BM _1 B T ) A" 

+ r; + ‘ + r; 


(4.5.5) 


where 

r" +1 = S 2 BM -1 Q n+1 + SB g] , (4.5.6) 

while the forward Euler formula 


A n+1 = A n + SX n+1 


gives 


(el + S 2 BM -1 B t ) A n+ ^ 


eA n + r" +1 


(4.5.7) 
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These update formulas are easily adopted into the generalized co- 
ordinate solvers discussed in the previous section. For each time step, in- 
tegration of the generized coordinates is first performed. The solutions axe 
used to update the right-hand side vector Q, which is then input into the 
Lagrange multiplier update module. The solution of the Lagrange multi- 
pliers is then used to correct the right-hand side vector Q with the current 
constraint forces B r A. The procedure can then be advanced to the next 
time step. 

4.6 Numerical Examples 

The computational techniques, namely the staggered multibody dy- 
namics solution procedure combining the generalized coordinate integrator 
presented in Section 4.3.2, the orientation update procedure presented in 
Section 4.4, and the constraint force solver presented in Section 4.5 with the 
computational procedure for the beam interned force presented in Chapter 
III, have been implemented into a Fortran 77 software package. The result 
is a robust method which solves the present formulation of the equations of 
motion of an arbitrary assemblage of flexible beams and rigid bodies. A dis- 
tinct feature of the present work is the computational preservation of total 
energy for undamped systems; this is obtained via the combined develop- 
ments of the objective strain increment/stress update procedure presented 
in Chapter III with the dynamic solution procedures of the present chap- 
ter. In order to demonstrate the current software capabilities, the following 
examples highlighting the flexible motion of the beam component are pre- 
sented. 

Example Spin-up maneuver. The first example is included to 
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verify the geometric stiffening phenomena exhibited by a rotating beam 
[15,23,30,33]. In this case, a pinned beam is subjected to the prescribed 
angular rotation 


6(t) = 


6 r / 2 15 2 . 2;rf ...... 

— [ 1 (cos 1) rad 0 < t < 15 sec 

15 L 2 2tt v 15 n 

( 6t — 45) rad t > 15 sec 


about the Z axis at the pinned end as shown in Figure 4.1a. The time 
history of the tip deflection relative to a reference frame coinciding with the 
prescribed angular position and the time history of several configurations 
of the beam are given in Figure 4.1b. An overall steady rotation of the 
beam gives rise to a centrifugal force which is responsible for a change in the 
bending stiffness that cannot be predicted using linear deformation theories. 
After initial increasing tip deflections, the beam begins to stiffen as the 
angular velocity increases due to the centrifugal inertia force. As the angular 
velocity reaches a constant state, the beam then reaches a steady state phase 
of small vibrations. The frequency of these vibrations is higher than those 
of a nonrotating beam due to the geometric stiffening effect. This example 
shows the capability of the nonlinear strain formulation to automatically 
account for this geometric stiffening, and the results are comparable to those 
presented in [33,100]. To reproduce these results with alternative methods 
as the substructuring technique [30], a convergence analysis based on the 
selection of mode shapes must be performed. 

Example ^.2 Planar motion of pinned beam. This next example 
demonstrates the combined large deformation and large rotation capabili- 
ties of the present formulation. Two different simulations are presented in 
which a pinned beam is subjected to the given initial velocity impulses shown 
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in Figures 4.2a and 4.2c. The resulting time histories of several deformed 
configurations due to the particular initial velocity impulse axe shown in 
Figures 4.2b and 4.2d. It is noted the versatility of the formulation in its 
ability to capture the response to a variety of situations in which different 
fundamental modes of the beam are excited. The approach avoids the dif- 
ficulty of tailoring the selection of modes shapes of the flexible components 
to the given problem at hand. The repeatability of the deformation shapes 
through time is due to the invariance of the internal force computations 
to the overall rigid motion. This property of computational objectivity is 
further illustrated in Figures 4.2b and 4.2d in which the time history of the 
strain, kinetic, and total energy is given over four revolutions of the motion. 
The nature of the time integration and internal force algorithms are such 
that the conservation of energy is retained computationally, as seen by the 
fact that the total energy remains constant for a number of revolutions. 

Example 4 ■ $ Closed-loop chain in free flight. This example, also 
analyzed in [100], demonstrates the capability of the present approach to 
model the dynamics of flexible multibody systems. A closed-loop chain con- 
sisting of four flexible links interconnected by hinges is subjected to a force 
and a torque at one of the hinges to induce an overall tumbling motion of 
the chain. The problem data is detailed in Figure 4.3a, and the sequence of 
the motion is given in Figure 4.3b. The results coincide with those presented 
in [100]. It is noted that the joint connection can easily be accounted for 
from a finite element assemblage which leaves the rotational degrees of free- 
dom free at the hinge location. The method was used to verify the results 
obtained using the Lagrange multiplier solver of Section 4.5. 

Example 4-4 Spatial motion of pinned beam. This example presents 
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the capability of the formulation to model the combined effect of spatial 
rotations with large deformations. As in Example 4.2, two different simu- 
lations are presented in which the pinned beam is subjected to an initial 
velocity impulse in the horizontal X-Y plane as shown in Figures 4.4a and 
4.4e. The beam is also subjected to a gravitational force acting along the 
-Z axis. The history of the spatial rotation and deformation throughout 
several revolutions axe shown in Figures 4.4b - 4.4c and 4.4e - 4.4f. The 
energy time histories for both simulations, as given in Figure 4.4g, again 
verify the computational objectivity of the algorithm for spatial motions. 

Example 1^.5. Spatial motion of double pendulum. This example 
presents the motion of a spatial double pendulum. The double pendulum 
is modeled with two beams; a spherical joint connects the last node of the 
first beam to the first node of the second beam and also pins the first node 
of the first beam. As shown in Figure 4.5a, the pendulum is subjected to a 
gravity field in the vertical Z-direction and an initial velocity impulse in the 
horizontal X-Y plane such that solely rigid motion is excited. The problem 
is run for four cases of increasing beam flexibility as given in Figure 4.5a. 
The spatial trajectories of the mass center of the second beam as projected 
on the X-Y and X-Z planes for each of the beam flexibilities are given in 
Figure 4.5b. The trajectory of the first case coincides exactly with a rigid 
body solution to the problem, and the slight deviation of the trajectories 
due to the increasing flexibility can be seen. The energy time histories for 
each case given in Figures 4.5c - 4.5f verify the computational conservation 
of energy. Again, the invariance of the internal force calculations in the 
three dimensional environment is witnessed by the negligible strain energy 
contribution within all of the flexible cases. The time integration of the 
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spatial kinematics preserves the balance between the kinetic and potential 
energies of the problem. 

It is noted that as the beam becomes more flexible, increasing 
amounts of strain energy slightly come into effect. Cases in which the beam 
was made even more flexible were run, and in these cases the conservation 
of energy was not maintained. After closer examination, it was determined 
the straining of the beam had reached unrealistic proportions beyond the 
elastic limit and the parameters were also unrealistic. An examination of 
the energy history provides an immediate detection of potentially invalid 
simulations. 

Example J^.6. Spatial motion of double pendulum. As a final ex- 
ample, the flexible double pendulum is given an initial velocity impulse to 
excite deformation motion as well as the rigid motion as shown in Figure 
4.6a. The first beam is given an impulse in the X-Y plane, while the second 
beam is given an impulse in the X-Z plane. Note the material parameters 
have been changed from the above example, and the gravity force has been 
removed. The resulting time histories of several deformed configurations axe 
given in Figure 4.6b, and the energy time histories are given in Figure 4.6c. 
The above example is included to show that the energy is still conserved for 
computations involving spatial deformations of the double pendulum. 

The inclusion of the energy history plots in some of the above ex- 
amples is a significant detail. This tool provides an effective method for 
analyzing the accuracy of the computational procedure. Numerous concep- 
tual and computational ‘bugs’ were detected from examination of the energy 
history. As such, it is believed this detail provides somewhat of an accu- 
racy measure of the given simulation. The figures referred to in the above 
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examples follow. 

4.7 Concluding Remarks 

A computational procedure suitable for the solution of equations of 
motion for flexible multibody systems has been presented. The equations 
of motion for an arbitrary assemblage of flexible beam and rigid body com- 
ponents can be constructed in a systematic manner as a result of the beam 
formulation presented in Chapter II. The dynamics axe effectively treated 
by the following developments in multibody dynamic solution procedures. 
A two-stage modification of the central difference algorithm is presented to 
integrate the translational coordinates and the angular velocity vector. This 
staggered from of the central difference method was adopted for multibody 
applications due to complications arising from the unavailability of the gen- 
eralized velocity vector at the time step at which the acceleration vector is 
evaluated. Given the solution to the angular velocity from the two-stage 
algorithm, the angular orientation is then obtained from the application 
of an implicit integration algorithm to the Euler parameter/angular veloc- 
ity kinematical relation. The constraint conditions, which are augmented to 
the formulation via Lagrange multipliers, axe enforced via a technique which 
implicitly integrates an alternative stabilized companion differential equa- 
tion for the Lagrange multipliers. The present multibody dynamics solution 
procedures axe effectively combined with the present internal force compu- 
tational procedure to achieve a computational preservation of total energy 
for undamped systems as demonstrated via several numerical examples. 

The next chapter presents the extension of the present methodology 
to model the dynamics of deployment /retrieval of the flexible members. 
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Figure 4.1a Spin-up maneuver: problem data. 
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Figure 4.2a Planar motion of pinned beam (1): problem data. 



xlO 5 Planar Motion of Pinned Beam: Four Revolution! 



Figure 4.2b Planar motion of pinned beam (l): motion history, energy history. 
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Figure 4.2c Planar motion of pinned beam (2): problem data. 




Figure 4. 2d Planar motion of pinned beam (2): motion history, energy history. 
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Figure 4.3a Closed-loop chain in free flight: problem data. 
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Figure 4.4a Spatial motion of pinned beam (l): problem data. 
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Figure 4.4d Spatial motion of pinned beam (2): problem data. 
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Figure 4.4e Spatial motion of pinned beam (2): motion history. 
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Figure 4.4g Spatial motion of pinned beam: energy history (1). 






xlO 5 Spatial Motion of Pinned Beam (2) 
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Figure 4.4h Spatial motion of pinned beam: energy history (2). 



nitial Beam Position vs. Initial Velocity Profile Material Properties 
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Figure 4.5b Spatial motion of double pendulum: projected trajectories. 
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Figure 4.5c Spatial motion of pinned beam: energy history (1). 



Spatial Double Pendulum: Case 2 
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Figure 4.5d Spatial motion of pinned beam: energy history (2). 




Spatial Double Pendulum: Case 3 
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Figure 4.5e Spatial motion of pinned beam: energy history (3). 




Spatial Double Pendulum: Case 4 
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Figure 4.5f Spatial motion of pinned beam: energy history (4). 
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Figure 4.6a Spatial motion of double pendulum: problem data. 
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Figure 4.6c Spatial motion of double pendulum: energy history. 






CHAPTER V 


DEPLOYMENT DYNAMICS OF BEAM STRUCTURES 
5.1 Introduction 

The purpose of this chapter is to use the present spatial beam formu- 
lation to model the dynamics of the beam’s deployment. Such a capability 
for simulating the extrusion of beam-like structures can offer a general util- 
ity to study various practical problems such as tethers, space assembly, hot 
rolling, and fluid jets. Much of the relevant literature addressing this type of 
problem has been motivated by space industry applications. Current space- 
craft designs employ flexible deployable appendages such as antennaes or 
solar arrays. The appendages, which are compactly stored during a launch 
phase, may be extended during the spacecraft’s attitude acquisition phase. 
As the relatively long appendages tend to be deployed with relatively small 
extension rates, the transient effect of deployment may be felt over a long 
period of time and is critical to the mission success. 

Initial works investigated the effect of appendage deployment on 
satellite control by modeling the appendages as point masses [101] and rigid 
bodies [102-103]. The deployment of flexible appendages from specific con- 
figurations of a spinning spacecraft was analyzed numerically in [104]. More 
general formulations for studying the effect of the deployment of beam-type 
appendages on the attitude stability of a satellite axe given in [105-10S], 
and a similar analysis for studying the deployment of plate-like structures is 
given in [109]. The above analyses employ appropriate modal coordinates to 
model the flexible deformations. As such, the classical vibrations approach 
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which employs a linear combination of spatial modal functions weighted by 
time-dependent generalized coordinates to model elastic deflections must be 
altered to account for the time- dependent spatial domain. As this approach 
may be somewhat questionable, an alternative method models the beam 
as a series of elastically connected rigid links, and then considers the links 
outside the rigid body at a given time in deriving equations of motion [110]. 
The equations of motion concerning the nonrotating dynamics of axially- 
moving strings [111-116], beams[117], and nonlinear elastica [118] have been 
formulated using an approach based on continuum mechanics. In the flavor 
of these last works, the present chapter analyzes the nonrotating dynamics 
of deployment. 

An axially moving beam-like structure is somewhat similar to a flow 
problem. The deployment may be thought of as a moving boundary problem 
in which the spatial domain occupied by the continuum is a function of time. 
Moving boundary problems are common in fluid dynamics, as typical heat 
conduction or diffusion problems involve phase changes from solid, liquid, 
or vapour states at an interface whose position is an unknown function of 
time. The fluid dynamics community has generated a great deal of research 
on numerical solution techniques for moving boundary problems which has 
provided much insight for the present work. 

A good account of numerical solution methods developed for moving 
boundary problems in heat flow or diffusion is given in [119-120]. In dis- 
cretizing the changing spatial domain encountered in these problems, either 
a fixed or a variable grid method may be adopted. The fixed grid method 
maintains a constant spatial grid sizing throughout the simulation. When 
the nodes are fixed in space, the location of the moving boundary does 
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not necessarily coincide with a node. Special numerical techniques must 
then be incorporated to locate the position and model the physics of the 
boundary [121-122]. Alternative approaches avoid these complications aris- 
ing from the unequal grid size near the moving boundary. In one method, 
the unequal interval is transferred to a more tractable region by moving 
the whole uniform grid system with the velocity of the moving boundary 
[123]. Another method varies the time step such that the boundary moves 
the distance of one fixed space mesh during that time [124-126]. Alter- 
natively, a deforming numerical grid which maintains an equal number of 
variable length space intervals may be employed to solve moving boundary 
problems [127]. This approach has a greater intuitive appeal as the moving 
boundary always corresponds to a specific nodal point. When employing 
this approach, the grid deformation must theoretically be accounted for in 
a proper manner. Within the context of the finite element method, the 
grid deformation can be accounted for by making the finite element basis 
functions implicit functions of time [128]. In this case, the time derivatives 
of the approximate displacements acquire additional convective- type terms 
within dynamic analyses. 

For dynamic analyses, space-time finite elements can be employed 
to formulate a time integration algorithm which accounts for the grid de- 
formation within a variable grid method. Previous works in the literature 
have used this approach with success to solve both parabolic and hyperbolic 
problems such as the one and two-dimensional heat equation of the Stefan 
ice/water interface problem [129-130], the Euler equations of compressible 
flow [131], and the advection-diffusion problem [132]. To solve flow prob- 
lems with time-varying boundaries in this manner, space-time variational 
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statements must first be defined. The simplest way to get the weak form 
of the governing equations is through the standard Galerkin finite element 
method [133-135]. However, the straightforward application of the Galerkin 
principle may result in a coupling of all the discrete steps in the time domain 
in the same manner that spatial nodal variables become coupled through 
a standard spatial finite element discretization. Obviously, this is highly 
undesirable for an efficient transient integration algorithm. To alleviate this 
problem, the discontinuous Galerkin method in time has been applied to 
space-time finite element formulations [136-138]. A more attractive method 
is to use Hamilton’s Law as the variational source for finite element dis- 
cretization procedures in the time domain [139-141]. In addition to leading 
to a step-by-step integration procedure, the use of a true variational law 
of mechanics as that of Hamiltonian as opposed to constructed principles 
as the general weighted residual methods provides a physical basis for the 
discretization procedure [139]. 

Given the above overview, the large deformation beam formulation 
presented in Chapter II is extended to model the axial deployment of the 
beam. The deployment is modeled by referring the motion variables to the 
changing spatial configuration of an undeformed beam extending with the 
prescribed deployment speed. Deforming finite elements are used to dis- 
cretize this changing spatial volume, and the internal force formulation of 
Chapter II and computation of Chapter III are retained. Due to the chang- 
ing spatial reference for the dynamic variables, the inertia operator acquires 
terms representing the convective rate of change of the variable in addition 
to those representing the local rate of change of the variable. To simplify 
the numerical treatment of the inertia operator, a Hamiltonian variational 
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formulation is introduced as the b&sis for a space-time finite element dis- 
cretization procedure. This space-time discretization of the Hamiltonian 
formulation results in a transient integration scheme which accounts for the 
deforming spatial grid. In this manner, an effective computational method 
has been developed to model the deployment of present beam formulation. 

The rest of the chapter will be organized as follows. Section 5.2 
presents the variable spatial grid description of the beam kinematics. The 
coupling between the convective effects and the temporal differentiation 
within the inertia operator are discussed in detail. The Hamiltonian for- 
mulation of the moving boundary value problem is derived in Section 5.3. 
The space-time discretization of the resulting Hamiltonian is carried out in 
Section 5.4, and care is exercised to annihilate unwanted spurious mecha- 
nisms in the discrete equations of motion. The computer implementation of 
the present approach and preliminary results on a planar inverse-spaghetti 
problems is illustrated in Section 5.5. 

5.2 Deployment Kinematics 

The problem under consideration, as shown below in Figure 5.1, is 
the axial extension of a beam-like structure from a stationary guide. As the 
beam extends from the guide, the spatial volume which deforms is constantly 
changing. To model this effect, the beam can be discretized by a growing 
number of finite elements of a fixed length or by an equal number of finite 
elements growing in length. These two approaches are illustrated in Figure 
5.2. It is seen that in the first approach, a special element of varying length 
must be used at either the boundary or the free end. This inconsistency 
is avoided in the second approach as the grid points are equally spaced. 
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** X 

Figure 5.1 Axial deployment subject to gravity. 

Conceptually, this variable or moving grid approach is more appealing, and 
it can readily be adopted into the present beam formulation. Theoretically 
and computationally, it remains to properly account for the effects of grid 
deformation. To this end, the kinematics of the beam deployment are de- 
scribed as follows. 

The present formulation is based on the use of an inertial reference 
for the beam kinematics. As detailed in Chapter II, the beam configuration 
is completely characterized using a position vector locating the neutral axis 
of the beam from the inertial origin and an additional reference frame which 
orients the cross-section from the inertial frame. To retain this concept and 
model the deployment of the beam, the displacement variables locating the 
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Figure 5.2 Fixed vs. Variable Grid Approaches 
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neutral axis are measured from a set of moving nodes which represent the 
position of a phantom beam rigidly moving from the guide in a prescribed 
manner. As shown in below in Figure 5.3, 



Figure 5.3 Deployment kinematics. 


the position vector r locating an arbitrary point within the changing spatial 
configuration occupied by the beam is described as 

r = { X(t) + u } T e + e T b . (5.2.1) 

In the above equation, X(t) represents time-dependent inertial coordinates 
of the phantom beam neutral axis. In reference to a particular moving node 
i, the position vector is given as 

r; = { Xi(t) + Ui } T e + £ t bi . (5.2.2) 

The location of the moving reference nodes X{(t) are determined by spa- 
tially dividing the length corresponding to a rigid deployment with an equal 
number of elements. If the prescribed deployment speed is a constant value 
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of V, the total length of the rigid phantom beam at time t is V t . The 
moving reference nodes are equally positioned between the free end and the 
boundary, and the position and velocity of the nodes can be prescribed as 

Xi (t) = ( - ~ 1 - V t (5.2.3) 

nett 

d -w 3 = {l ^r v (5 - 2 ' 4) 

where title represents the fixed number of deforming elements. The kine- 
matic relation (5.2.2) is thus an Eulerian description as material points of 
the beam are not specifically tracked. 

As the displacement variables are measured from a moving reference, 
the material time derivative definition must be introduced to obtain the 
velocities and accelerations of a material particle. The definition of the 
material time derivative is given as [45] 


D_ _ d_ d 
Dt ~ dt Vj dxj 


(5.2.5) 


where x j are the inertial coordinates of a spatial point within the deformed 
beam configuration, and Vj axe the instantaneous velocity coordinates at 
that point. For use within the present beam formulation, this definition 
must be transformed such that the spatial derivative is taken with respect 
to the convected coordinates This convected coordinate definition is 
given as 


D 

Dt 


d „ d 

s ft + ViTii w i 


(5.2.6) 


where the rotational tensor Tij maps the inertial coordinates to the con- 
vected coordinates. It is recalled that the convected reference frame was 
introduced in Chapter II to obtain a more physically effective deformation 
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description. As the internal force is computed from the convected spatial 
coordinates, the transformation (5.2.6) is introduced to achieve a consistent 
computation of the inertia force. For the particular Eulerian description 
(5.2.1), the appropriate material time derivative becomes 


® .l + yi 

Dt at d( 

V* — T\\ V\ + T \ 2 V 2 + T13 U3 


(5.2.7) 

(5.2.8) 


as the translational variables and the cross-section orientation are solely 
functions of the neutral-axis length coordinate £. To complete this defini- 
tion, an expression for the instantaneous velocity u* is obtained from the 
kinematic definition (5.2.2). From this definition, the spatial coordinates 
are given as 

X{ = X{ + Ui + Rji £j (5.2.9) 

where R^ is the rotational tensor corresponding to the cross-section orien- 
tation. The instantaneous velocity of this point becomes 


Vi — Vi -|- ili -|- Rji VJjfc 


Vi 


dXj 

dt 



(5.2.10) 

(5.2.11) 


From this expression, the velocity V* of (5.2.8) can be written as 


V* = u* + 0 * T £ 


(5.2.12) 


where 


ii* — T\j ( Vj + itj ) 

* T 

0* = {0 , 5 *i 1 UJ3 — 531 u?i , 5 2 1 ^1 ~ ^2 } 
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correspond to the translational and rotational contributions respectively as 
derived from (5.2.8) and (5.2.10). 

The velocity and acceleration of a material particle can now be 
derived by applying (5.2.6) to the Eulerian description (5.2.1). For rotating 
coordinate systems, the definitions of the time and space derivatives for the 
body fixed components of a given vector are given as 


d e _ d b 
dt dt 

dd ~ di 

The velocity of a material 


+ w , u> — - 

+ K , K — — 

particle thus becomes 


— r / 
dt 

(5.2.13) 

dR T 

~3( R ' 

(5.2.14) 


D r 
~Dt 


(V + u ) + V* 


d e u 

~w 


+ 


u 


+ V* k 


(5.2.15) 


and likewise the acceleration 


D 2 


Dt 2 


r = e T 


< + 2 V* + V’ £? + 


dt 2 


dt dd 


dd 2 


d e n dv i ai- 

' dt ; 


dd 

r d b u 

~df 


+ UlCb + 


dd 


+ 


* , d b Cj d b k „ _ _ _ x 

v ( "ar + ar +u, ' c + '“ j) + 

v ( + KK ) + 

dV* „„ dV* . 

-5T k + y 


(5.2.16) 


The equations of motion as derived from the stated beam kinematic descrip- 
tion will be discussed next. 
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5.3 Hamiltonian Equations of Motion 

The beam formulation presented in Chapter II is based on the well- 
known Cauchy equilibrium equations of motion for a solid continuum. The 
principal of virtual work, which is a spatial variational form of the partial 
differential equations of motion, is given as 


6F 1 + 6F S = 6F e + 6F t (5.3.1) 


where the inertia operator 6F 1 , the internal force operator SF s , the exter- 
nal force operator 8F E , and the traction operator 6F T are identified from 
(2.3.1). Due to the beam kinematic assumptions, the traction operator is 
identically equal to zero and will henceforth be neglected. This variational 
form provides a basis for the implementation of spatial finite element ap- 
proximation methods to lead the partial differential equations of motion to 
a set of ordinary differential equations in time. Standard numerical integra- 
tion techniques can then be employed to obtain a solution to the equations 
at discrete time levels. 

For the deployment problem, a variational form of the temporal op- 
erators is considered such that a finite element discretization in time can be 
applied simultaneously with the discretization in space. As the principle of 
virtual work is used as the basis for the spatial discretization, Hamilton’s 
Law is used as a variational source for this space-time finite element dis- 
cretization. Hamilton’s law is deduced from an integration of the principle 
of virtual work over an arbitrary time interval given as [98] 



8F 1 dt + 



/: 


6F S dt 


SF e dt . 


(5.3.2) 
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An integration by parts is then performed on the inertia term 

r *2 r l 2 

SF 1 dt = 

'ti 

to lead the equation (5.3.2) to the form 


r 8F 1 di = r f P fi 8n 

Ju Jti Jo 


dd dt 


(5.3.3) 


f I { p r, Sri - Oij feij ~ 6 r ifi } dd dt = 
Jtr J# 

[ J p riSri < 


*2 
J tl 


(5.3.4) 


A well-known form of the above variational statement is obtained by requir- 
ing that the arbitrary virtual changes 6r{ vanish at the time limits t\ and 
t 2 as 

Sri ( ti ) = Sri ( t 2 ) = 0 


With this assumption, the boundary term on the right-hand side of (5.3.4) 
can then be neglected such that 



(ST + SW)dt 


0 


where T represents the kinetic energy and SW represents the virtual work of 
the applied and internal forces. This interpretation is known as Hamilton’s 
principle. In general though, one can treat r,( t\ ) and r%( t 2 ) as undeter- 
mined quantities. The distinction between Hamilton’s Law and Hamilton’s 
principle is thus whether these temporal boundary terms are retained or 
neglected. 

This distinction between Hamilton’s Law and Hamilton’s principle 
becomes important in the numerical approximation of the variational state- 
ments. Although the boundary term is often seen as irrelevant in deriving 
the equations of motion, it is important for numerical solution techniques. 
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The boundary term must be retained to achieve correct approximate solu- 
tions from a Hamiltonian variational statement. Early works which imple- 
mented a space-time discretization of Hamilton’s principle resulted in an 
inconsistent treatment of initial data [133-135,142]. In contrast, the initial 
conditions can be properly incorporated into the variational formulation of 
Hamilton’s Law. Additionally, when Hamilton’s Law is used as the varia- 
tional statement, Co continuous approximating functions may be employed 
such that an attractive step-by-step integration procedure results [139-140]. 

To interpret the well-known Hamilton’s Law for the present deploy- 
ment problem, an integration by parts must be performed on the inertia 
term 

[ SF 1 dt = f f p Sri dd dt (5.3.5) 

Jt x Jti Jd(t) Dt 

which contains the material time derivative definition (5.2.6). The proper 
interchange of the time derivative and a time dependent spatial volume 
integral is given by the Reynolds transport theorem as [45] 


»L fF * - /*>'*-* • 

where p is- the mass density and F is an arbitrary function of the material. 
From this theorem, an integration by parts can be performed on (5.3.5) to 
yield 




D r, 
Dt 


Sri dd 


D ri D Sri 
~Dt DT 


<2 


dt) dt 


(5.3.6) 


The interpretation of Hamilton’s Law for Eulerian continuum formulations 
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is thus given as 



D Sri 
Dt 


<7 ij S£ij 



Srifi } dd dt = 


(5.3.7) 


The above variational statement provides an attractive source for a 
concurrent space-time discretization procedure for the present deployment 
problem. The complex acceleration expression (5.2.16) need not be explicitly 
considered as the temporal weak form requires only the velocity expressions. 
In addition, consistency with the internal force operator is achieved as the 
spatial derivatives contained within the inertia axe reduced to first order. 
Finally, the deforming spatial grid which has been introduced to model the 
changing spatial volume occupied by the beam will be taken into account 
automatically. 

To complete the formulation, an expression for the inertia operator 
is derived by incorporating the kinematic expressions into (5.3.6); the details 
of the derivation are given in Appendix A. By neglecting coupling terms 

* T 

caused by the rotational contribution 0* £ within the velocity V* given in 

(5.2.12), the interior term of the inertia operator becomes 


sf '- = r i 

Jtl Jv 


D Vi D 8r l 
Dt Dt 


dd dt 

du '\ T 

f dSu 

. .* dSu \ 

+ ti* 



+ “ aT ) 

6u + 

u * Sk 

| dt dt . 

(5.3.8) 


In this expression, 


m 


= [ pdA , J = f peF dA 

JA Ja 
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represents the beam mass per unit length and cross-sectional inertia prop- 
erties, and 


Cj 6a 
k 6a . 

represent the angular velocity and curvature variations, respectively. Like- 
wise, the resulting expression for the boundary term of the inertia operator 
becomes 


6ui 

6k 


d 6a 

~6T 

d 6a 

~w 


+ 


+ 



These expressions (5.3.8) and (5.3.9) are much more tractable to handle 
numerically than the acceleration term given earlier in (5.2.16). 

To complete (5.3.7), the nonlinear internal force operator, as derived 
in Chapter II, is given as 

[’ 6F S = r j <*“ T iaT ) I B ] T {^} d(dt , (5.3.10) 
and likewise the external force operator is given as 




d£ dt 


(5.3.11) 


5.4 Space-Time Discretization 

In conventional finite element models, the displacement field u is 
approximated with a linear combination of spatial interpolation functions 
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and time-dependent nodal displacements as given by (3.3.1). This concept 
can be modified to account for mesh deformation by making the spatial in- 
terpolation an implicit function of time. Thus when node motion is allowed, 
the appropriate spatial approximation is given as 

npe 

u( u) =s £ N, ( m ) u, (() . (5.4.1) 

1=1 

For a simultaneous discretization in time, the nodal displacements are ap- 
proximated in the same manner as a linear combination of temporal inter- 
polation functions and nodal values at discrete time steps as 

npe 

uj(t) ~ Y Nj ( t ) u/ J . (5.4.2) 

j= l 

From (5.4.1) and (5.4.2), the space-time discretization is given as 

npe npe 

u ~ u = Ni (0 Nj ( t ) u/ J . (5.4.3) 

l=i J=i 

In this expression, the value u/ J corresponds to the value of u at the spatial 
node I and time tj. The nodal value is thus the displacement with respect 
to the node Xj(tj) whose position is prescribed throughout the simulation 
as a function of the deployment speed. 

5.4.1 Discrete Epilations of Motion 

The space- time discretizations of the inertia, internal force, and ex- 
ternal force operators within Hamilton’s Law (5.3.7) are given as follows. 
The following analysis is specialized to the case of planar motion where the 
appropriate degrees of freedom become 


u 


{ tii,U2,0 } 
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corresponding to two translational displacements which locate the neutral 
axis and a single angle of rotation which orients the body-fixed reference 
frame. These three degrees of freedom can be treated in the same manner. 
The inertia expression reduces to the form 


6Fj 

SF b 


f = r j 

Jt i 

t - { j * 


D 6u T ,, Du 

M 


Dt 


Dt 


8u 


Dt 


M 8u d£ 


(5.4.4) 

(5.4.5) 


where 


M = Diag { m , m , J } 

This expression is evaluated by substituting the approximations 

+ 

d Nj (0 


Dt 
D 8u 


/= i j=i 

npc npc 


U 


tj 


Dt 




Nj ( 0 


d£ 


6uj 


tj 


1=1 7=1 

into (5.4.4) and (5.4.5). The result for the interior term is written as 


SFi 1 = { 6u tl 8u ts } 


An A12 

A 2 i A22 


u 


u 


(5.4.6) 


where the spatial element kernels constituting the A matrices are given in 
Appendix B. Likewise, the boundary terms of the inertia operator, evaluated 
as 


8F, 


' = / 
J«.h) 


D u(t 2 ) 


Dt 


Su(t 2 ) d£ 




DJjt 1 ) 

Dt 


8u(t \ ) d£ 


is approximated using 


^ (o + u «•*■* 


J =1 

npc 


8u(t u t 2 ) = 


J= 1 


+ u 
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where it is reiterated that 


u = ( ui,u 2 ,6 ) 

« *'* = T u ( V\ + ) + T 12 ( V 2 + u 2 t2 ) . 


The discrete form of (5.4.5) is given as 
SFf, 1 = { 6u fl 6u h } 


Bi 

0 

C, 

0 


0 

B 2 

0 

c 2 


{£} 

{:*) 


+ 


(5.4.7) 


where the spatial kernels B\, B 2 ,C\, and C 2 are also given in Appendix B. 
The discrete form of the linear internal force operator is given as 



6F S dt 


{ 6u u 6u h } ^ 


K 

K 



(5.4.8) 


where K is the linear Timoshenko beam stiffness matrix and h is the time 
step between t\ and t 2 . This linear discretization is presented to illustrate 
the general flow of the algorithm; the extension of the algorithm incorpo- 
rating the nonlinear internal force is discussed in Section 5.4.3. Finally, the 
external force is discretized as 


J‘‘ SF e dt = { 6u" 6u" + f,[ } ( 5 - 4 - 9 > 

In all the above discretizations, the double integration over the 
changing spatial domain and the arbitrary time interval has been performed 
using reduced integration methods. The reduced integration methods were 
first chosen as a similar application implementing a space-time discretization 
of Hamilton’s Law concluded that the reduced integration methods resulted 
in unconditional stability of the transient integration algorithm. In con- 
trast, the full integration methods resulted in a conditional stability of the 
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transient integration algorithm [141]. A discrete Fourier analysis, presented 
in Section 5.4.2, concludes that the reduced integration gives a consistent 
approximation of the original partial differential equations of motion. 

From the discretizations (5.4.6 - 5.4.9), the transient integration 
algorithm is given as follows. The discrete equations of motion arising from 
the space-time procedure are written as 




+ C x u‘‘ 
+ C2 u* 7 



- f(/‘ l 

- f(/‘ l 


+ f t2 )\ 

+ )J 


(5.4.10) 

(5.4.11) 


for the linear system. Given the initial conditions u tl and ti* 1 


the first block 


of the above equations can be solved for u <J as 


{ ^12 - ~ K } u h = - { A u . (5.4.12) 

It can be shown from the explicit expression of A X2 given in Appendix B 
that the left-hand side matrix in the above equation is positive definite. 
Thus from a solution to (5.4.12) for u‘ 2 , the second block equation can then 
be solved for it 1 7 as 


B 2 u t7 = { A 2X - j A' } + { A 22 - ^ K - C 2 } u h -f ~( /*> + f t7 ) . 

Again, it can be seen that B 2 is also positive definite. These solutions u l 2 
and u t2 form the initial conditions for the next time step, and thus the time 
finite elements result in a step-by-step integration formula. 
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5.4.2 Consistency Analysis 

A discrete Fourier analysis provides a mechanism for examining the 
consistency of the discrete equations of motion with respect to the origi- 
nal partial differential equations [143-144]. The strong form of the partial 
differential equations of motion are given as 


+ 2V^ + = 0 . 


(5.4.13) 


dt 2 ' dxdt dx 2 

As the main desire is to analyze the characteristics of the inertia operator, 
the following approximations have been made in the above equation. The 
convective velocity V* has been approximated by the deployment speed V , 
and only axial vibrations have been retained such that V 2 = V 2 4- c 2 
where c 2 represents the elastic wave speed E / p. 

The traditional Fourier analysis seeks a general harmonic wave so- 
lution to the above equation of the form 


u 


Un e i{wi - kx ) 


(5.4.14) 


where w is the circular frequency, k the wave number, and i the imaginary 
number i = y/—l- The substitution of (5.4.14) into (5.4.13) yields the 

characteristic equation 

{ - w 2 + 2V kw - V 2 k 2 } u 0 = 0 . (5.4.15) 

The fundamental relation between the frequency and the wave number as 
dictated by this characteristic equation is thus given by 

w * = ( V- c* ) k * 


(5.4.16) 



where non-dimensional parameters 
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w* = wh , k* = kh , V* = -j- , c * = -j 

have been introduced for an arbitrary step-size h and element length £, 
respectively. 

For the discrete equations of motion, the difference equations are 
obtained by assembling the spatial nodes of two interior beam elements 
with two sequential steps in time. The spatial nodes, designated as j — 
1 > j , 3 + 1, are spaced a distance £ apart, and the time nodes, designated 
as i — 1 , i , i + 1, are spaced a time step of h apart such that it*' — ► ttj , . 
For the finite element equations resulting from the reduced integration of 
the variational equations, the difference equations become 

£ r 

[ U i+1,»+1 - 2 Uj+1,; + Uj + i,i_i 4- 2 Uj'i+1 - 4 Uj t i + 

2 - 2 uj — l,t + ] + 

2 V r , 

~ 7 ~ [ U i+l,i+l - «> — l.i-t-l - J + 

47 t u i+i,«+i - 2 -f- 4- 2 Uj+i,j - 4 4- 

2 Uj-i,i 4 — 2 ] 

= 0 . (5.4.17) 


Likewise, the difference equations for a fully integrated version become 


6 h 
2 V 


4 



[ «>+i,H-i — 4 Uj + i'i + Uj+i'i-i + 2 uj ii+ i — 8 uj'i 4- 

2 Uj.j-i + - 4u/'-l,i + U;-!,,-! ] + 

[ u i+i,i+i — Uj-i.i+i — 4- ] -f 

[ u j+ i,«+i — 4 Uj,i+i 4- + 2 Uj + 1,, — 8 uji + 
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2 u j-l,i + u j+l,i-l ~ 4 Uj,i - 1 + ttj-l.f-l ] 


= 0 


(5.4.18) 


To compare the characteristics of the difference equations to the 
characteristics of the continuum equation given in (5.4.15), a harmonic so- 
lution to the difference equations of the form 


u 




e i(wt-kz) 


is sought. By substituting the harmonic solution into the difference equa- 
tions (5.4.17) and (5.4.18), the characteristic equations of the discrete cases 
become 


{i 


£ 2 ~ ^ Tr s'mk£ sin wh 

- 2V — ~h~ 

+ V 2 { 1 - ^ w 2 } k 2 

1 4 J 


0 (5.4.19) 


for the reduced integration version (5.4.17) and 


- 2 r 2 , , h 2 _ 2 1 l 2 

1 3 V 4 3 h? V 2 


2 V 


sin kl sin wh 


} 


w 


4 f 1 V* ! 

+ 3 { h 2 + 02 ^ 


i 2 


(5.4.20) 


for the fully integrated version (5.4.18) where 


k = 


sinf 

£ 

2 


sin 


wt 


W = 


t£ 

2 


for both cases. 

It is now observed from a careful comparison of the characteris- 
tic equation for the continuum case (5.4.15) to those for the discrete cases 
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(5.4.19) and (5.4.20), the latter corresponding to the full integration scheme 
results in an inconsistent representation. Spurious mechanisms emanate 
from nonphysical rigid-body motions. The possibility of introducing spuri- 
ous mechanism by the preceding finite element discretization can be deter- 
mined from the solution of the above equations with w = 0. The reduced 
integration characteristic equation (5.4.19) reduces to 

k 1 2 = 0 =*> k = n = 0,1,2, (5.4.21) 


for this analysis. As the highest admissible deformation mode shape admit- 
ted by linear shape functions is k / £, the only potential solution of (5.4.19) 
is k — 0. This solution coincides with the correct rigid-body mode solution 
of the continuum characteristic equation (5.4.16). The latter case of full 
integration leads to 


1 - k 2 - i — — 1 

3 h 2 V 2 + 3 Wi 2 + £ 2 * 


0 


which does not admit a physically valid solution. It is also noted that the 
representation (5.4.20) does not converge in the limit to the true character- 
istic equation (5.4.15), whereas the representation (5.4.19) does have this 
property. Finally, the dispersion relation of the reduced integration differ- 
ence equations can be derived as 


w* -4- k* 

tan y = ( V* - c* ) tan — 

to be compared to the true dispersion relation (5.4.16). This is done below in 
Figure 5.4 where it is seen the approximation is valid for the range k* < j. 
The reduced integration version is thus the consistent discretization of the 
variational equations of motion. 
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Figure 5.4 Dispersion relations. 


5.4.3 Nonlinear Solution Strategy 

To include the nonlinear force term within the procedure, the equa- 
tions 


- f U B r- M * - it >} 

are solved via the Newton-Rhapson iteration technique. The displacement 


An 
A 21 


A 12 
Ai2 


is obtained first by iteratively solving a nonlinear version of (5.4.10). To 
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this end, a solution to the following residual equation 

r = A u u* 1 + A 12 u * 2 - ['* [ [ B ] T a dd dt - f rh * 1 

Jti Jv 

= 0 (5.4.22) 

is obtained by iterating about the linearized set of equations 


E(k+i) Au* 2 (fc+1) — -r(* +1 ) . (5.4.23) 


The above equations are solved at the ( k + l) tA iteration for incremental 
displacements Au* 3 which update the solutions obtained at the k th iteration 
via 

u (/t+i) = u (k ) 3 + Au * 2 . 

The iterative procedure is started with = u* 1 and continued until 

a convergence criterion is reached. The solution matrix E introduced in 
(5.4.23) is obtained from (5.4.22) as 


E = A 12 - £ ( K g + Km ) 

where K G and K M correspond to the material and geometric stiffness ma- 
trices of the nonlinear internal force. 

To complete the procedure, an evaluation of the time integral of 
the internal force within the residual (5.4.22) is necessary. Using a reduced 
integration interpretation, this is given as 



8F S dt 


{ 6u ** 8u* 3 



+ a' 2 ) ) \ 
+ a* 2 ) ) j 


where the internal force is evaluated using an average of the displacements 
between the two time steps. 
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It is important to note that in the moving node formulation the 
displacement coordinates u" and are defined with respect to different 

spatial grids; the coordinates it" are defined from the nodes X n whereas 
the coordinates u" +1 are defined from the nodes X" +1 . Therefore, prior 
to the computation of the displacement average for the evaluation of the 
internal force, the coordinates u n must be extrapolated to the grid X n+1 . 
In addition, the internal force computation discussed in Chapter III was 
based on an incremental procedure in which the current stress state was 
obtained by updating a past stress state with an increment of stress. This 
stress increment is a function of the displacement increment A u between a 
current configuration u and the past configuration u n as 

Au = u n+1 - it" . 

To incorporate this concept within the moving node formulation, the dis- 
placement increment must be also be computed from variables defined on a 
consistent grid. 

To this end, an extrapolation of the variables u" to the grid X n+l 
must be defined in a manner consistent with the present internal force for- 
mulation. The stress update remains unaffected by the extrapolation of 
it" if the constant elemental strain states of the past configuration are rep- 
resented without alteration. A consistent extrapolation thus must retain 
the orientation of the elemental convected reference frames. By retaining 
the convected reference frame orientation as determined from the original 
translational coordinates x n for the extrapolated displacements x n ' , the 
deformed position of the beam neutral axis can equivalently be described 
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with either coordinates as 


x n = X n + u n 
x n = X n+1 + u nl . 

Figure 5.5 below illuminates the above concept. For a given element 
bounded by nodes i and i + 1, the aj axis of the convected reference frame 
is determined from 


*£j+1 Xi 

e 


l = II Xi + 1 - Xi I . 


I 


f 



Figure 5.5 Extrapolation of displacement coordinates. 

To retain the orientation of the element convected frame, the extrapolations 
must be such that 

•j'j+l •J'i+l X{ 


e> 


i 


(5.4.24) 
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Given a constant deployment velocity V , the nodal positions X, are auto- 
matically prescribed such that the length between two nodes 

= Xi + 1 — X, = — 

nele 

is a function of time as given by (5.2.3), It is easily seen that, starting from 
the first node which remains fixed at the boundary of the guide as 

V'n v n + 1 V f 

Aj — A j — A j , 

the rest of the nodes for i — 2, nele + 1 can be generated as 

xUi = x; + c x 

= X! + ( Jf i+ . - * ) • 

If the extrapolated displacements are generated from previous values in the 
same sequential manner for nodes i = 2, nele + 1 as 

u[ — Ul 

= Uj + ) 5 

it is easily shown the condition (5.4.24) is satisfied. The position of the 
deformed neutral axis and thus the strain state is unaltered by the extrap- 
olation. The rotational variables are extrapolated in the same manner such 
that a constant curvature state is retained without alteration. These extrap- 
olations are performed on variables prior to an internal stiffness computation 
and are otherwise not a part of the integration algorithm itself. 
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5.5 Results 

The computational techniques, namely the space- time finite element 
integration algorithm combined with the internal force computation dis- 
cussed in Chapter III, have been implemented into a Fortran 77 software 
package. The software was first verified by repeating planar examples re- 
ported in Chapter IV using a fixed node reference. The formulation was 
then tested on an analysis of a sheet of paper issuing from a rigid horizontal 
guide into a uniform gravitational field and compared to results reported 
in [118]. In work of [118], the nonlinear equations of motion of an elastica 
that moves out of a horizontal guide at a constant velocity are shown to 
depend on two parameters, namely a dimensionless weight- to-stiffness ratio 
H = mg£* / El and a dimensionless velocity v = V£y/ml / ~EI. In these 
expressions, m is the mass per unit area of the paper, g is the acceleration 
of gravity, l is the length of the paper, El the bending stiffness, and V the 
constant velocity of the paper ejection. The analysis was performed with 
H = 50 and v 2 = 100,50,20,10 corresponding to V = 92,65,41,20 
in/sec respectively. The trajectories of the beam tip for the various deploy- 
ment speeds are shown in Figure 5.6 and compared to the static shape of a 
fully extended cantilevered beam with given properties such that fi = 50. 
The results are comparable to those reported in [18]. Figure 5.7 gives an 
animated history of the beam deformations throughout various stages of 
the deployment. The deformation shapes are again compared to the static 
shape. 






Figure 5.7 Motion history for various deployment speeds (V). 
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5.6 Concluding Remarks 

An effective computational method has been developed for the mod- 
eling of beam deployment which retains the advantages of the present beam 
formulation. A deforming spatial grid has been introduced to model the 
changing spatial volume. Due to the changing spatial reference for the 
dynamic variables, the inertia operator acquires terms representing the con- 
vective rate of change of the variable in addition to those representing the 
local rate of change of the variable. A variational formulation is presented 
in Hamiltonian form to simplify the numerical treatment of the inertia op- 
erator. A space-time discretization of the Hamiltonian formulation is im- 
plemented which accurately accounts for the deforming spatial grid in the 
transient integration scheme. Reduced integration techniques were used to 
analyze the required integrations over both the space and time domains; 
a discrete Fourier analysis concluded that this discretization resulted in a 
consistent approximation to the governing partial differential equations of 
motion. To effect the moving node formulation within the present internal 
force computational procedure, the solutions of a past configuration are ap- 
propriately extrapolated to the grid representing the current configuration. 
The results obtained from the present computational procedure agree with 
those reported in [118]. 

The next chapter summarizes the work presented in this disserta- 


tion. 



CHAPTER VI 


CONCLUSIONS 


6.1 Summary of Work 

The present thesis work has focused on the dynamic analysis of 
three-dimensional elastic beams which experience large rotational and large 
deformational motions. A realistic mathematical model of a spatial flexible 
beam was developed as an integral kernel of a general multibody dynam- 
ics methodology. The model accounts for both large rotations and large 
deformations as typically experienced by a flexible component within an ar- 
ticulated structure. Computational solution techniques were then derived 
and implemented which enhance and exploit various features inherent in the 
formulation. The resulting methodology provides a tool to study the effects 
of component flexibility on the performance of multibody systems. 

To model the spatial dynamics of highly flexible beams, the beam 
motion is described using an inertial reference for the translational displace- 
ments and a body-fixed reference for the rotational quantities. Finite strain 
rod theories are then defined in conjunction with the beam kinematic de- 
scription. The resulting strain measures account for the effects of stretching, 
bending, torsion, and transverse shear deformations, and thus model poten- 
tial deformations of typically lightweight and highly flexible space structure 
appendages. In addition, due to the kinematic description of the beam mo- 
tion, the form of the equations of motion are similar to that of rigid body 
dynamics. As such, these equations can be solved with numerical solution 
procedures developed for general multibody dynamic systems. Another ad- 
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vantage of the present formalism is the use of a convected coordinate rep- 
resentation of the Cauchy stress tensor and a conjugate strain definition to 
model the beam deformation. This is in contrast to many finite-deformation 
analysis which typically adopt the Piola-Kirchoff stress tensor to model the 
deformation. As such, the stress/strain representation of the present formu- 
lation coincides with the actual strains measured by sensors rotating with 
and operating on the deformed structure. An easy interface is thus possi- 
ble with control-type applications such as vibration suppression and slewing 
maneuvers. 

The numerical treatment of the beam formalism is considered in de- 
tail in the present thesis. A computational procedure for the internal force 
has been derived from the continuum formulation. As the formulation is 
based on an inertial reference for the beam dynamics, the degrees of free- 
dom of the flexible component contain information of both rigid and flexible 
deformation motions without distinction. The rigid motions must not in- 
fluence the strain computation. To this end, a unique computation which 
exploits characteristics of the strain formulation as well as finite rotation 
theory was derived to filter out the rigid body motions embedded within 
the degrees of freedom. The procedure was proven to remain invariant to 
arbitrary finite rigid rotations of the beam while accurately modeling the 
beam strain. 

The present internal force computation was successfully interfaced 
with multibody dynamic solution procedures. As a consequence of the 
present beam formulation, the structure of the equations of motion is iden- 
tical for both rigid and flexible components of an arbitrary multibody con- 
figuration. The following numerical solution procedures, which have been 
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developed and successfully implemented in the present work, axe thus appli- 
cable for both the rigid and flexible components of the articulated structure. 
An interlaced application of the central difference algorithm is used to inte- 
grate the translational coordinates and the angular velocity vector. Given 
the solution for the angular velocity, an implicit algorithm is used to dis- 
cretize the Euler parameter/angular velocity kinematical relation to obtain 
the corresponding angular orientation. The constraint conditions, which 
are augmented to the formulation via Lagrange multipliers, are enforced via 
a separate procedure which implicitly integrates an alternative stabilized 
companion differential equation for the Lagrange multipliers. The numeri- 
cal techniques were originally tested on multirigid body systems and then 
successfully extended to models incorporating the present beam formula- 
tion. The combined developments of the objective internal force compu- 
tation with the dynamic solution procedures result in the computational 
preservation of total energy for undamped systems. This distinct feature 
has been demonstrated in several examples. 

The final development presented in this thesis is the modeling of 
the dynamics of deployment/retrieval of the beam formulation. To model 
the deployment, a moving spatial grid is employed as a reference for the 
dynamic variables. This reference corresponds to the configuration of a 
deployed beam as if it were rigid. The introduction of this moving spatial 
reference leads to complexities in the beam inertia operator. For this reason, 
a Hamiltonian variational formulation is employed to simplify the numeri- 
cal treatment of the inertia operator. Space-time finite elements are used 
to discretize the Hamiltonian formulation, resulting in a transient integra- 
tion scheme which accurately accounts for the deforming spatial grid. The 
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methodology is successfully interfaced with the internal force computational 
formalism, thus retaining the large rotation/large deformation modeling ca- 
pabilities of the present work. 

The major contributions of this work axe summarized as follows: 

a) A formulation for the spatial motion of flexible beams has been de- 
veloped based on physically applicable strain definitions represen- 
tative of those measured from actual sensors located and operating 
on a rotating flexible structure. 

b) A numerical procedure for the accurate computation of the inter- 
ned force has been developed from the continuum formulation and 
proven to be invariant to arbitrary rigid motions of the beam. 

c) Multibody dynamic solution procedures have been developed and 
applied to models of articulated structures incorporating the present 
beam formulation. 

d) The dynamics of deployment/retrieval have been modeled by incor- 
porating a moving node reference for the present beam formulation 
and a space-time discretization of the corresponding Hamiltonian 
variational statement. 

e) The computational procedures developed to simulate the dynam- 
ics of multibody systems with flexible beam components have been 
implemented into a software testbed, providing a mechanism for 
dynamic analysis of space structures, robot assemblies, and similar 
applications. 

6.2 Directions for Further Research 


Accurate numerical simulations of flexible multibody systems can 
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be utilized to facilitate interdisciplinary research. The techniques used in 
designing control schemes for robotic type applications can be enhanced 
with this simulation capability. As the implementation of a maneuvering 
strategy developed from rigid models is not guaranteed to be successful, it 
is important to determine to what extent the flexibility affects the robot 
performance. 

A potential application for this type of simulation capability is the 
in-space construction facility being developed for the construction, main- 
tenance, and repair of future massive space structures. The main body 
of the crane is designed to carry the structural components to be assem- 
bled. A mobile transporter which controls an attached remote manipulator 
system moves along the crane arm to provide fine positioning and delicate 
readjustment tasks required of the assembly process. Limited numerical ex- 
periments which have been performed on flexible models of the space crane 
suggest that the development of control strategies and maneuvering speeds 
need be further addressed, and the effect of flexibility on the tip positioning 
accuracy need be determined. Control strategies that work when the sys- 
tem configurations change in a manner effecting rapidly varying frequency 
responses must be developed for the successful operation of the space crane. 

To effect on-board adaptive control strategies, real-time software 
simulation capabilities must also be developed. To this end, parallel algo- 
rithms that address the highly nonlinear equilibrium equations augmented 
with complex kinematic and control constraints must be derived and imple- 
mented within current multi-processor technologies. 

In addition to applications involving robotics and control, another 
area of research is the prediction of conditions which lead the present beam 
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model to unstable motion. The influence of the deployment of flexible ap- 
pendages on the spacecraft attitude stability may be parameterized in terms 
of the deployment velocity and the orbital angular velocity. In addition, it 
should be explored whether the present computational procedure can pre- 
dict the onset of chaotic motion within the nonlinear beam model due to 
a change in certain input conditions. In general, the present flexible beam 
computational methodology provides a tool for further analysis of many 
engineering operations. 
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APPENDIX A 


It is desired to obtain an expression for the variational operator 
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corresponding to the inertial operator of the Hamiltonian formulation. The 
velocity of a material particle was given in (5.2.15) as 
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where V* is the neutral-axis or £ component of the instantaneous velocity 
of a spatial point as defined in (5.2.8). The virtual displacements Sr are 
defined as a set of infinitesimally small displacements of purely kinematic 
nature that are consistent with the system constraints at a given instant 
of time. Therefore, for the deployment problem Sr is defined for a specific 
instant of the moving reference where X is held fixed and the definition 
(2.2.10c) 
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remains valid. The following expression is thus obtained 
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by applying the material time derivative operator (5.2.7) including the rotat- 
ing coordinate transformations (5.2.13 - 5.2.14) to the virtual displacement 
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definition. With this information, (A.l) can now be evaluated as 
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where the well known expressions 
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have been introduced. The above is integrated throughout constant cross- 
sectional area coordinates represented by l. Typically, the last two terms 
of the above expression which axe linear in i vanish as an integration is 
performed over symmetrical cross-sections. In this case however, V* is also 
a function of t 

V* = u* + 0* T t 


as derived in (5.2.12), and additional terms are generated. The final expres- 
sion is given as 
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represents the beam mass per unit length and cross-sectional inertia prop- 
erties and 

r = { o , 0 2 * j 2 , 0 3 * h } T . 

In the same manner, the boundary term is evaluated as 


^ = [j 


— L 

{ ( V + u* ) + tit* Tjj } T m 8u + { u> T + u* k t } J 8a + 

^ R t 8at + r T k t R 8 u d£ 

9( 


It is seen that the convective velocity V* gives rise to a slight coupling 
between the translational and rotational degrees of freedom. For the present 

T 

analysis, we will neglect the effect of 0* t within the convective velocity 

V*. 



APPENDIX B 


The space-time discretization of the inertia operator within Hamil- 
ton’s Law was given in (5.4.4) as 
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for the interior term and 
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was given in (5.4.7) for the boundary term. In the above terms, the A, B , 
and C matrices have been derived using both full and reduced integration 
techniques in analyzing the integrations of the changing spatial domain and 
the arbitrary time interval inherent in the formulation. The results axe given 
below in terms of individual spatial finite element kernels to be assembled 
across the spatial domain via the manner standard of finite elements. 

The A matrices axe given explicitly as 
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The spatial element kernels as derived from reduced integration techniques 
are given as 
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where 


( ti * + u 2 * ) 


is the average nodal velocity u * for a given element. The same kernels 
integrated exactly are given as 


Mi 

p A l 

' 2 

1 ' 


— * 





6 

_ 1 

2 


M 2 s 

P A 

r — ! 

i ? ', * _ 

* * 

- 1“ 

• • 


4 



il2_ 




3 “ 

3 

m 2 a 

II 

Ts 

lu 

tol 

a 

' o - 

l 

-1 ' 
0 

m 3 

P A 

t 2 

l 

-1 

= — 7- u 

a 

1 

-1 




Finally, the factors 





1 

4 

1 

4 


correspond to the exact vs. reduced integration analysis respectively. With 
this information, the difference equations (5.4.14) and (5.4.15) to be ana- 
lyzed with the discrete Fourier procedure are easily derived. 


Likewise, the spatial kernels of the boundary term are defined as 


B i 
B 2 


= B ( ) , C x = C ( ti* ) , C 2 = C ( Tij Vj ) 

'Bit*') + A u*' Tu D A Tu D 0 

= A u 2 t2 Tn D B ( €*» ) + Au 2 Ij T 12 D 0 

A 6 *' Tu D A 9 *' T 12 D B ( t *' ) 
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where 


A Ui 


Uil - Uii 


i 


and the matrices are defined as 


B = 
C = 
D = 


pAi r i i 

4 l 1 1 



1 

1 


for the reduced integration evaluation and 


B 

C 


p A t 
6 

pA 

6 


D = 


1 

6 



2iiJ + uj 
u\ + 2 u* 2 


for the exact integration evaluation. 





